STEVEN SCOX 


MATRIX 
ANALYSTS 


Matrix Analysis 


Collection Editor: 
Steven J. Cox 


Table of Contents 


Preface 
1.1. ‘Preface: tio, Matrix Analysis® 2020.4 -9h0eeacekoees teen cd ei eg cae eda eA Ped eee 2 
Matrix Methods for Electrical Systems 
2.1 Nerve Fibers and the Strang Quartet ....... 0.0... cece cece eee ee cee ene ee ees 5 
2.2 CAAM 335 Chapter 1 Exercises ............0..20 020 c cece eee en eect ete e eens 15 
SOlUItIONS. ja tecidectere's salen seein aerated Seg gael eee ae ap bake eee ee eee ate P oe et ee 17 
Matrix Methods for Mechanical Systems 
Sel: “Al Uniaxial Truss.) ce8 ie aeteiad ote nee heh eee e he ee tele lense tele Rta 19 
3:2 AvSmall Planar-Truss:.2424.00 3.40.4 fa cect aca odie eases een tad oie vacide ere eee ede eae 24 
3.3 The General Planar Truss ............... 0002 cece eee eee eee eee ene teen eee e eee eeeees 27 
3.4 CAAM 335 Chapter 2 Exercises ............0..20 0c eee e eee eee eee ne dente eeeees 30 
The Fundamental Subspaces 
AL .- Column: Space 22 sk Choe Ph te ks BEE ONG ge hee Boe ee BE Se Ueto Seed BEES 33 
AD RNa ES Pace? iin Bre oka Bae Pata PR Bh a eee Ate OLS Ns Sh ent Cienega deen ee eG hie tedsh ets 34 
4.3. The Null and Column Spaces: An Example ................2002 0 ccc cece cece eee ee eee eeteeeees 36 
AA Nett NulliSpace skcc 2. et aoe as bg tae Macys Moa eg ete eee Oe Bl ye ea ee nila ee We AO 
ALS) ROW Space: .o2 tt ee ce eho gta pheitaaehe eid Pena ctd gata eee a aa eho 2 eee de 41 
4.6 Exercises: Columns and Null Spaces ............. 0.0.2 c cece eee nee eect tenet teen eens 42 
4.7 Appendices/Supplements ............0.0 cece cece eee e cette eee eee e eee nent eeeees 42 
Least Squares 
Bid“ Lie@astsSquares: «02 ty cent kee kind aay okay eelee ab thie inane oh ca aed oe ae eee ee A7 
Matrix Methods for Dynamical Systems 
6.1 Nerve Fibers and the Dynamic Strang Quartet ............0..2 0002 e eee cece cee eee eee 53 
6.2° The Laplace Transform’ +: :2..0¢.cea¢oceeide ca ceate 4 eee ed eee a va ete eed ea ceedasaeee 58 
6.3. The Inverse Laplace Transform ................00 cece eee eee een e eee e eee e teenies 60 
6.4 The Backward-Euler Method ............... 0.000 cece e eee ene n eee cee eee ee eees 61 
6.5 Exercises: Matrix Methods for Dynamical Systems ................2.0-2 2c e eee cee eee eee 62 
6:6\-Supplemental':2.3.3.53323.053 8 Dl ae hs oe Olt a eee baed Meb ane abhl Labakaece, 63 
Complex Analysis 1 
7.1 Complex Numbers, Vectors and Matrices .............0.. 2c cece cece eee eect ete e ee eees 75 
(22° ‘Complex*Functions 3.2. 30 eee eke pee ee ee ae ee ee 77 
723. Complex Ditterentiation .j...245f.05 ois el ees eg gee Se ee oe 79 
7.4 Exercises: Complex Numbers, Vectors, and Functions .............2.. 0200s eee eee eee eee ees 84 
Solutions. 24) 3.2 e's gee he sees poets Satie des edi awe lane ey eu ee eee ene eee 85 
Complex Analysis 2 
8:1. “Cauchy’s Theorem 9.02 ged eid tend. teed ated ote oe ee a nate bait cate chee eka 87 
8.2 Cauchy’s Integral Formula ............. 0 eee eee eee cence tence eens eens eneeeaeee 89 
8.3 The Inverse Laplace Transform: Complex Integration ................0.02 200 eee een eee ees 93 
8.4 Exercises: Complex Integration ......... 0.0... cece ec c nc n cnet eee e eet e een eens 94 


The Eigenvalue Problem 


9.1. Introduction: 2.124.220. he eR ee es es cd 95 
9.2. “Phe: Resol vents 2220s. aur gcils ke dig onh geal oe aol oxen pent ees weeea ce eee none Bedi aneeth els 96 
9.3 The Partial Fraction Expansion of the Resolvent ................ 0.20 cee eee cence ee eee eee eees 97 


9.4 The Spectral Representation ............. 200. e eee eee tee tenet nett e eee eeeees 99 


9.5 The Eigenvalue Problem: Examples ............... 00. c cece cece n eee eee nee ee eee ees 100 

9.6 The Eigenvalue Problem: Exercises ............... 0. c cece eee eee cnet nent eee eee nees 102 
10 The Symmetric Eigenvalue Problem 

10.1 The Spectral Representation of a Symmetric Matrix ............0.. 0002 eee eee eee eens 103 

10.2 Gram-Schmidt Orthogonalization ......... 0.0... ence ene nee teen teen ees 107 

10.3 The Diagonalization of a Symmetric Matrix .......... 0.0. c cee cece eee eee eens 108 
11 The Matrix Exponential 

DEST SOVverviews 256222 ches cette ehh a ae eh et eh ee 111 

11.2 The Matrix Exponential as a Limit of Powers ................ 2.0. e ccc eee eee tenet eees 113 

11.3. The Matrix Exponential as a Sum of Powers ................ 2200 e eee e eee e ee tenet n eens 114 

11.4 The Matrix Exponential via the Laplace Transform ..................2.200 2c cee eee eee eeees 116 

11.5 The Matrix Exponential via Eigenvalues and Eigenvectors .................0. 02 eee eee eee 117 

11.6 The Mass-Spring-Damper System .............0 0.2 ccc e eee eee eect eee e eee nent eeees 121 
12 Singular Value Decomposition 

12.1 The Singular Value Decomposition .......... 00.0... e cece cece nnn cence cence eens 125 
Glossary 0556 cis ses nee ee ees BE ee ba eon Se ee Paes Meld Rte Wee ed 130 
Bibliography: \ioe.ess eet oeecens ee oes Pee kere ee ee ee ee ee Eo Se eee eat ae 134 
TnG OX? ys se easel dete Meh ite Rhee See ea ei eee eee eaten laine tage ina ee eee at eaten? 136 


2 CHAPTER 1. PREFACE 


Chapter 1 


Preface 


1.1 Preface to Matrix Analysis’ 


Matrix Analysis 


Figure 1.1 


!This content is available online at <http://cnx.org/content /m10144/2.8/>. 


Bellman has called matrix theory ’the arithmetic of higher mathematics.’ Under the influence of Bellman 
and Kalman, engineers and scientists have found in matrix theory a language for representing and analyzing 
multivariable systems. Our goal in these notes is to demonstrate the role of matrices in the modeling of 
physical systems and the power of matrix theory in the analysis and synthesis of such systems. 

Beginning with modeling of structures in static equilibrium we focus on the linear nature of the relation- 
ship between relevant state variables and express these relationships as simple matrix-vector products. For 
example, the voltage drops across the resistors in a network are linear combinations of the potentials at each 
end of each resistor. Similarly, the current through each resistor is assumed to be a linear function of the 
voltage drop across it. And, finally, at equilibrium, a linear combination (in minus out) of the currents must 
vanish at every node in the network. In short, the vector of currents is a linear transformation of the vector 
of voltage drops which is itself a linear transformation of the vector of potentials. A linear transformation 
of n numbers into m numbers is accomplished by multiplying the vector of n numbers by an m-by- n ma- 
trix. Once we have learned to spot the ubiquitous matrix-vector product we move on to the analysis of the 
resulting linear systems of equations. We accomplish this by stretching your knowledge of three-dimensional 
space. That is, we ask what does it mean that the m-by- n matrix X transforms ¥” (real n-dimensional 
space) into ’? We shall visualize this transformation by splitting both ” and R™ each into two smaller 
spaces between which the given X behaves in very manageable ways. An understanding of this splitting of 
the ambient spaces into the so called four fundamental subspaces of X permits one to answer virtually 
every question that may arise in the study of structures in static equilibrium. 

In the second half of the notes we argue that matrix methods are equally effective in the modeling and 
analysis of dynamical systems. Although our modeling methodology adapts easily to dynamical problems 
we shall see, with respect to analysis, that rather than splitting the ambient spaces we shall be better 
served by splitting X itself. The process is analogous to decomposing a complicated signal into a sum of 
simple harmonics oscillating at the natural frequencies of the structure under investigation. For we shall see 
that (most) matrices may be written as weighted sums of matrices of very special type. The weights are 
eigenvalues, or natural frequencies, of the matrix while the component matrices are projections composed 
from simple products of eigenvectors. Our approach to the eigendecomposition of matrices requires a brief 
exposure to the beautiful field of Complex Variables. This foray has the added benefit of permitting us a 
more careful study of the Laplace Transform, another fundamental tool in the study of dynamical systems. 

—Steve Cox 
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Chapter 2 


Matrix Methods for Electrical Systems 


2.1 Nerve Fibers and the Strang Quartet’ 


2.1.1 Nerve Fibers and the Strang Quartet 


We wish to confirm, by example, the prefatory claim that matrix algebra is a useful means of organizing 
(stating and solving) multivariable problems. In our first such example we investigate the response of a nerve 
fiber to a constant current stimulus. Ideally, a nerve fiber is simply a cylinder of radius a and length / that 
conducts electricity both along its length and across its lateral membrane. Though we shall, in subsequent 
chapters, delve more deeply into the biophysics, here, in our first outing, we shall stick to its purely resistive 
properties. The latter are expressed via two quantities: 


1. p;, the resistivity in Qcm of the cytoplasm that fills the cell, and 
2. Pm, the resistivity in Qcm? of the cell’s lateral membrane. 


!This content is available online at <http://cnx.org/content /m10145/2.7/>. 


CHAPTER 2. MATRIX METHODS FOR ELECTRICAL SYSTEMS 


A 3 compartment model of a nerve cell 


Figure 2.1 


Although current surely varies from point to point along the fiber it is hoped that these variations are 
regular enough to be captured by a multicompartment model. By that we mean that we choose a number 


N and divide the fiber into N segments each of length +. Denoting a segment’s 


Definition 2.1: axial resistance 


R, = 28 
Ta 
and 
Definition 2.2: membrane resistance 
Rao 
2nay 


we arrive at the lumped circuit model of Figure 2.1 (A 3 compartment model of a nerve cell). For a fiber 
in culture we may assume a constant extracellular potential, e.g., zero. We accomplish this by connecting 


and grounding the extracellular nodes, see Figure 2.2 (A rudimentary circuit model). 


A rudimentary circuit model 


R. R; R 


i i 


Figure 2.2 


Figure 2.2 (A rudimentary circuit model) also incorporates the exogenous disturbance, a current 
stimulus between ground and the left end of the fiber. Our immediate goal is to compute the resulting 
currents through each resistor and the potential at each of the nodes. Our long-range goal is to provide 
a modeling methodology that can be used across the engineering and science disciplines. As an aid to 
computing the desired quantities we give them names. With respect to Figure 2.3 (The fully dressed circuit 
model), we label the vector of potentials 


ee ( V1 LQ %«®3 ) 
and the vector of currents 
y= ( Yr Y2 %Y3 Ya Ys Ye I 


We have also (arbitrarily) assigned directions to the currents as a graphical aid in the consistent application 
of the basic circuit laws. 
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The fully dressed circuit model 


Figure 2.3 


We incorporate the circuit laws in a modeling methodology that takes the form of a Strang Quartet|1]: 


(S1) Express the voltage drops via e = — (Az). 

(S2) Express Ohm’s Law via y = Ge. 

(S3) Express Kirchhoff’s Current Law via A’y = —f. 
(S4) Combine the above into A7GAz = f. 


The A in (S1) is the node-edge adjacency matrix — it encodes the network’s connectivity. The G 
in (S2) is the diagonal matrix of edge conductances — it encodes the physics of the network. The f in (S3) 
is the vector of current sources — it encodes the network’s stimuli. The culminating A7G‘A in (S4) is the 
symmetric matrix whose inverse, when applied to f, reveals the vector of potentials, x. In order to make 
these ideas our own we must work many, many examples. 

2.1.2 Example 
2.1.2.1 Strang Quartet, Step 1 


With respect to the circuit of Figure 2.3 (The fully dressed circuit model), in accordance with step (S1) (list, 
1st bullet, p. 8), we express the six potential differences (always tail minus head) 


€1 = %1 — &2 
€2 = 1% 
€3 = 12 — @ 


€4 = V3 


€5 = %3— %4 


€6 = &4 
Such long, tedious lists cry out for matrix representation, to wit e = — (Ax) where 
-1 1 0 0 
0 -1 0 O 
0 -1 1 0 
A — 
0 oO -1 O 
0 oO -1 1 
0 0 oO -1 


2.1.2.2 Strang Quartet, Step 2 


Step (S2) (list, 2nd bullet, p. 8), Ohm’s Law, states: 
Law 2.1: Ohm’s Law 
The current along an edge is equal to the potential drop across the edge divided by the resistance 
of the edge. 


In our case, 


yy = Ges F=1,3,5 and y= Zr, 7 =2,4,6 


or, in matrix notation, y = Ge where 


SOO fo Oo Fh 
ooo fo + ° 
Soo OFF CO OO 
o° P+ ee) 
oe COO Oo 
+ OS: SO 2S: 3S 


2.1.2.3 Strang Quartet, Step 3 


Step (S3) (list, 3rd bullet, p. 8), Kirchhoff’s Current Law”, states: 


Law 2.2: Kirchhoff’s Current Law 
The sum of the currents into each node must be zero. 


In our case 
10 = Yi = 0 


Yi — y2— y3 = 0 


¥3 — Y4— Y5 = 0 
2"Kirchhoff’s Laws": Section Kirchhoff’s Current Law <http://cnx.org/content /m0015/latest /#¢current> 
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Y5 — Yo = 0 
or, in matrix terms 
By=—f 
where 
—-1 0 0 0 0 to 
1 -1 -1 0 
B= and f= 
0 0 —-1 -1 


2.1.2.4 Strang Quartet, Step 4 
Looking back at A: 


-1 1 0 60 
0 -1 0 O 
ja: 0 -1 1 0 
0 O -1 0O 
0 oO -1 1 
0 O O -Il 


On substitution of the first two into the third we arrive, in accordance with (S4) (list, 4th bullet, p. 8), at 
A’ GAg = f. (2.1) 


This is a system of four equations for the 4 unknown potentials, x, through x14. As you know, the system 
(2.1) may have either 1, 0, or infinitely many solutions, depending on f and A’GA. We shall devote (FIX 
ME CNXN TO CHAPTER 3 AND 4) to an unraveling of the previous sentence. For now, we cross our 
fingers and ‘solve’ by invoking the Matlab program, fibl.m ° . 


Shttp://www.caam.rice.edu/~+caam335/cox/lectures/fibl.m 
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Results of a 64 compartment simulation 


Potential along 64 compartment fiber 


"1 T — T T 
10F 4 
ob 4 
— & 
=> 
E 
* of | 
8 
SE “f 
4 1 4 n 1 4 n 4 
0 0.1 o2 o3 u4 os o6 o7 o8 og 


z (cm) 


Figure 2.4 


Results of a 64 compartment simulation 


Axial current along 64 compartment fiber Membrane current along 64 compartment fiber 


y (BA). 
py 


> yo 0t8F 4 


2 (em) ~~ z (cm) 


Figure 2.5 


This program is a bit more ambitious than the above in that it allows us to specify the number of 
compartments and that rather than just spewing the x and y values it plots them as a function of distance 
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along the fiber. We note that, as expected, everything tapers off with distance from the source and that the 
axial current is significantly greater than the membrane, or leakage, current. 


2.1.3 Example 


We have seen in the previous example (Section 2.1.2: Example) how a current source may produce a potential 
difference across a cell’s membrane. We note that, even in the absence of electrical stimuli, there is always a 
difference in potential between the inside and outside of a living cell. In fact, this difference is the biologist’s 
definition of ‘living.’ Life is maintained by the fact that the cell’s interior is rich in potassium ions, K*, 
and poor in sodium ions, Na‘, while in the exterior medium it is just the opposite. These concentration 
differences beget potential differences under the guise of the Nernst potentials: 


Definition 2.3: Nernst potentials 


T, |N T, |k 
— R [Na] and Ex = Flore 


where R is the gas constant, T’ is temperature, and F' is the Faraday constant. 


Associated with these potentials are membrane resistances 


Pm,Na and Pm K 


that together produce the p,, above (list, item 2, p. 5) via 


and produce the aforementioned rest. potential 


+ 


B= (— Ex ) 
m = Pm 
Pm,Na Pm,Na 


With respect to our old circuit model (Figure 2.3: The fully dressed circuit model), each compartment 
now sports a battery in series with its membrane resistance, as shown in Figure 2.6 (Circuit model with 
resting potentials). 
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Circuit model with resting potentials 


R. R, R 


i i 


Figure 2.6 


Revisiting steps (S1-4) (p. 8) we note that in (S1) the even numbered voltage drops are now 


€2 = £2 — Em 


€e4 = 43 — En, 


eg = 14— Em 
We accommodate such things by generalizing (S1) (list, Ist bullet, p. 8) to: 
e (S1’) Express the voltage drops as e = b — Ax where b is the vector of batteries. 
No changes are necessary for (S2) and (S3). The final step now reads, 


e (S4’) Combine (S1’) (p. 13), (S2) (list, 2nd bullet, p. 8), and (S3) (list, 3rd bullet, p. 8) to produce 
A?GAz = ATGb+ f. 


Returning to Figure 2.6 (Circuit model with resting potentials), we note that 


m 


and A’Gb= —” 


m 


Ee oO fF CO KF O&O 
a & 
Ke FP FF CO 
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This requires only minor changes to our old code. The new program is called fib2.m 4 and results of its use 
are indicated in the next two figures. 


Results of a 64 compartment simulation with batteries 


Potential along 64 compartment fiber 


~Hi + 


~~ 
= 
E 
** as} 
we 
ao 
6 1 1 1 
0 0.1 02 of of OS O08 OF 08 of 
z (cm) 
Figure 2.7 


“http://www.caam.rice.edu/~caam335/cox/lectures/fib2.m 
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Results of a 64 compartment simulation with batteries 


Axial current along 64 compartment fiber Membrane current along 64 compartment fiber 


< < 0.02 
0 1 4 2 (om) ? 0 dl 0 4 2 (cm) 
(a) (b) 
Figure 2.8 
2.2 CAAM 335 Chapter 1 Exercises’ 
Exercise 2.1 (Solution on p. 17.) 


In order to refresh your matrix-vector multiply skills please calculate, by hand, the product ATGA 


in the 3 compartment case and write out the 4 equations in the vector equation we arrived at in 
step (S4) (2.1): A7GAz = f. 


Exercise 2.2 

We began our discussion with the ’hope’ that a multicompartment model could indeed adequately 
capture the fiber’s true potential and current profiles. In order to check this one should run fib1.m® 
with increasing values of N until one can no longer detect changes in the computed potentials. 


e (a) Please run fibl.m’ with N = 8, 16, 32, and 64. Plot all of the potentials on the same 
(use hold) graph, using different line types for each. (You may wish to alter fib1.m so that 
it accepts N as an argument). 


Let us now interpret this convergence. The main observation is that the difference equation, (2.2), 
approaches a differential equation. We can see this by noting that 


acts as a spatial ’step’ size and that x,, the potential at (k — 1) d(z), is approximately the value of 
the true potential at (k — 1) d(z). In a slight abuse of notation, we denote the latter 


a ((k — 1) d(z)) 


>This content is available online at <http://cnx.org/content /m10299/2.8/>. 
Shttp://www.caam.rice.edu/~+caam335/cox/lectures/fibl.m 
Thttp://www.caam.rice.edu/~caam335/cox/lectures/fib1l.m 
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Applying these conventions to (2.2) and recalling the definitions of R; and R,, we see (2.2) become 


ma? —x (0) + 2x (d(z)) — x (2d (z)) i 2nad (z) 


x (d(z)) =0 
Pi d(z) Boe a 
or, after multiplying through by xray: 


Pi d (2?) 


. We note that a similar equation holds at each node (save the ends) and that as N — oo and 
therefore d(z) — 0 we arrive at 


x(z)=0 (2.3) 


e (b) With » = 2% show that 


x(z) =asinh (/2u2) + Bcosh (V2u2) (2.4) 


satisfies (2.3) regardless of a and (3. 
We shall determine a and @ by paying attention to the ends of the fiber. At the near end we find 


ma? x (0) — x (d(z)) 


a Od 
which, as d(z) — 0 becomes 
d pito 
ace =— 2. 
a Oa (2.5) 
At the far end, we interpret the condition that no axial current may leave the last node to mean 
d 
eis = 2: 
Sr () =0 (2.6) 


e (c) Substitute (2.4) into (2.5) and (2.6) and solve for a and @ and write out the final x (z). 
e (d) Substitute into x the 1, a, p;, and pm values used in fibl.m® , plot the resulting function 
(using, e.g., ezplot) and compare this to the plot achieved in part (a) (p. 15). 


Shttp://www.caam.rice.edu /~caam335/cox/lectures/fibl.m 


Solutions to Exercises in Chapter 2 


Solution to Exercise 2.1 (p. 15) 


2.1 Feedback 


The second equation should read 
=24- > 222 — @3 is v2 


R; Rm 


17 


(2.7) 
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CHAPTER 2. MATRIX METHODS FOR ELECTRICAL SYSTEMS 


Chapter 3 


Matrix Methods for Mechanical Systems 


3.1 A Uniaxial Truss’ 


3.1.1 Introduction 


We now investigate the mechanical prospection of tissue, an application extending techniques developed in 
the electrical analysis of a nerve cell (Section 2.1). In this application, one applies traction to the edges of a 
square sample of planar tissue and seeks to identify, from measurement of the resulting deformation, regions 
of increased ‘hardness’ or ‘stiffness.’ For a sketch of the associated apparatus, visit the Biaxial Test site ? . 


3.1.2 A Uniaxial Truss 


A uniaxial truss 


Figure 3.1 


As a precursor to the biaxial problem (Section 3.2) let us first consider the uniaxial case. We connect 3 
masses with four springs between two immobile walls, apply forces at the masses, and measure the associated 
displacement. More precisely, we suppose that a horizontal force, f;, is applied to each mj, and produces a 
displacement x;, with the sign convention that rightward means positive. The bars at the ends of the figure 
indicate rigid supports incapable of movement. The k; denote the respective spring stiffnesses. The analog 


'This content is available online at <http://cnx.org/content /m10146/2.6/>. 
*http://health.upenn.edu /orl/research/bioengineering/proj2b.jpg 
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of potential difference (see the electrical model (Section 2.1.2.1: Strang Quartet, Step 1)) is here elongation. 
If e; denotes the elongation of the jth spring then naturally, 


€y =X 


€2 = 1%Q2— @ 


€3 = %3— Xr 


€4 = —2X3 
or, in matrix terms, e = Ax, where 
1 O 0 
-1 
A — 
0 -1 1 
0 O -1 


We note that e; is positive when the spring is stretched and negative when compressed. This observation, 
Hooke’s Law, is the analog of Ohm’s Law in the electrical model (Law 2.1, Ohm’s Law, p. 9). 
Definition 3.1: Hooke’s Law 
1. The restoring force in a spring is proportional to its elongation. We call the constant of propor- 
tionality the stiffness, k;, of the spring, and denote the restoring force by y;. 
2. The mathematical expression of this statement is: y; = kje,;, or, 
3. in matrix terms: y = Ke where 


ky 0 0 0 

O. B00 
Kee : 

Or. Feo 

0 0 “Ok 


The analog of Kirchhoff’s Current Law (Section 2.1.2.3: Strang Quartet, Step 3) is here typically called 
‘force balance.’ 


Definition 3.2: force balance 
1. Equilibrium is synonymous with the fact that the net force acting on each mass must vanish. 


2. In symbols, 
yi — yo— fr =9 


yo — 3 — fo =0 


y3 — ya — fg = 0 


3. or, in matrix terms, By = f where 
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As in the electrical example (Section 2.1.2.4: Strang Quartet, Step 4) we recognize in B the transpose of 
A. Gathering our three important steps: 


e = Ax (3.1) 
y = Ke (3.2) 
ATy = f (3.3) 


we arrive, via direct substitution, at an equation for 7. Namely, 
(ATy a f) > (A? Ke = f) > (ATK Ax — 7) 


Assembling A’ K Ax we arrive at the final system: 


ky + ko —ko 0 XY fi 
—kp  kotkgs —kg rw |=] fe (3.4) 
0 —kg kg + ka x3 fs 


3.1.3 Gaussian Elimination and the Uniaxial Truss 


Although Matlab solves systems like the one above with ease our aim here is to develop a deeper understand- 
ing of Gaussian Elimination and so we proceed by hand. This aim is motivated by a number of important 
considerations. First, not all linear systems have solutions and even those that do do not necessarily possess 
unique solutions. A careful look at Gaussian Elimination will provide the general framework for not only 
classifying those systems that possess unique solutions but also for providing detailed diagnoses of those 
defective systems that lack solutions or possess too many. 

In Gaussian Elimination one first uses linear combinations of preceding rows to eliminate nonzeros below 
the main diagonal and then solves the resulting triangular system via back-substitution. To firm up our 
understanding let us take up the case where each k; = 1 and so (3.4) takes the form 


2 -1 0 Ly fi 
-1 2 -1|/| » |=| fh (3.5) 
0 -1l 2 x3 fs 


We eliminate the (2,1) (row 2, column 1) element by implementing 


1 
new row 2 = old row 2+ so 1 


bringing 
2 -l1 0O Ly fi 
0 3 —l v2 = fe + 4 
0 -1 2 X3 fs 


We eliminate the current (3,2) element by implementing 


2 
new row 3 = old row 3+ sow 2: 
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bringing the upper-triangular system 


2 -1 O Ly fi 
0 3 —1 v2 = fe + A 
00 4 x3 fat 7B+4 


One now simply reads off 
fit 2fo + 3fs 
L3 = fo 


This in turn permits the solution of the second equation 


2(mtht#) pAsropeh 
3 7 2 


6 


and, in turn, 
es tat fi _ 3fit2fe + fa 
2 4 
One must say that Gaussian Elimination has succeeded here. For, regardless of the actual elements of f, we 
have produced an x for which A? K Ax = f. 


3.1.4 Alternate Paths to a Solution 


Although Gaussian Elimination remains the most efficient means for solving systems of the form Sx = f it 
pays, at times, to consider alternate means. At the algebraic level, suppose that there exists a matrix that 
‘undoes’ multiplication by S in the sense that multiplication by 2~! undoes multiplication by 2. The matrix 
analog of 2-!2 = 1 is 

ses =T 


where J denotes the identity matrix (all zeros except the ones on the diagonal). We refer to S~! as: 


Definition 3.3: Inverse of S 
Also dubbed "S inverse" for short, the value of this matrix stems from watching what happens 
when it is applied to each side of Sx = f. Namely, 


(Se = f) = (S-1Se¢ = Sf) = (Im = Sf) = (ex = Sf) 


Hence, to solve Sx = f for x it suffices to multiply f by the inverse of S. 


3.1.5 Gauss-Jordan Method: Computing the Inverse of a Matrix 


Let us now consider how one goes about computing S~!. In general this takes a little more than twice the 
work of Gaussian Elimination, for we interpret 


Si tg=] 


as n (the size of S') applications of Gaussian elimination, with f running through n columns of the identity 
matrix. The bundling of these n applications into one is known as the Gauss-Jordan method. Let us 
demonstrate it on the S appearing in (3.5). We first augment S with I. 


S2 eal, ge a ae 
AU. Ai, LOA 
Oe Br coro A 


23 


We then eliminate down, being careful to address each of the three f vectors. This produces 


a ha D0 
0 # -1 | $ 10 

4 1 2 
De OS ge il 2g. 43d 


Now, rather than simple back-substitution we instead eliminate up. Eliminating first the (2,3) element we 
find 


2-10 {41 00 

09 0133 3 

00 ¢/|43 31 
Now, eliminating the (1,2) element we achieve 

2-0 0 [8 1 4 

Og Ol a oe 

00% |: 21 


In the final step we scale each row in order that the matrix on the left takes on the form of the identity. This 
requires that we multiply row 1 by $, row 2 by 3. and row 3 by 3, with the result 


100 |32 3 4 
O42, aed 
001,44 4 23 


Now in this transformation of S into I we have, ipso facto, transformed I to S~!; i.e., the matrix that 
appears on the right after applying the method of Gauss-Jordan is the inverse of the matrix that began on 


the left. In this case, 


= 


BIR NIF BI 
NRF Re Nile 
BIW NIF BIR 


One should check that S~!f indeed coincides with the 2 computed above. 


3.1.6 Invertibility 


Not all matrices possess inverses: 
Definition 3.4: singular matrix 
A matrix that does not have an inverse. 
Example 
A simple example is: 


Alternately, there are 


Definition 3.5: Invertible, or Nonsingular Matrices 
Matrices that do have an inverse. 
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Example 
The matrix S that we just studied is invertible. Another simple example is 


3.2 A Small Planar Truss’ 


We return once again to the biaxial testing problem, introduced in the uniaxial truss module (Section 3.1.1: 
Introduction). It turns out that singular matrices are typical in the biaxial testing problem. As our initial 
step into the world of such planar structures let us consider the simple truss in the figure of a simple swing 
(Figure 3.2: A simple swing). 


A simple swing 


Figure 3.2 


We denote by x; and x2 the respective horizontal and vertical displacements of m, (positive is right and 
down). Similarly, f; and fy will denote the associated components of force. The corresponding displacements 
and forces at mz will be denoted by x3, x4 and fs, f4. In computing the elongations of the three springs we 
shall make reference to their unstretched lengths, L,, Le, and L3. 

Now, if spring 1 connects {0,—L1} to {0,1} when at rest and {0,—L1} to {#1, 72} when stretched then 
its elongation is simply 


e1 = 2m? + (ma + a)’ — hh (3.6) 


’This content is available online at <http://cnx.org/content /m10147/2.6/>. 
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The price one pays for moving to higher dimensions is that lengths are now expressed in terms of square 
roots. The upshot is that the elongations are not linear combinations of the end displacements as they were 
in the uniaxial case (Section 3.1.2: A Uniaxial Truss). If we presume, however, that the loads and stiffnesses 
are matched in the sense that the displacements are small compared with the original lengths, then we may 
effectively ignore the nonlinear contribution in (3.6). In order to make this precise we need only recall the 

Rule 3.1: Taylor development of the square root of (1 + t) 
The Taylor development of 21+ t about t = 0 is 


Va +t=1+5+0(2) 


where the latter term signifies the remainder. 


With regard to e; this allows 


ey = 2212 t Lo? t 22511 t Eye Ty (3 7) 
= Tyyf21 4 aati? 4 Bee Ty 
v1 t+a27 x17 +a” Qu 
a = Ty + ie + 242+ 1,0 ( es + 2) Ty 
= a1? +a? L,O x47 +a2? 2x2 = 8) 
= @2 2D, 1 Tyee Tage 
If we now assume that 
vy? + x92 
rae is small compared to x2 (3.9) 


then, as the O term is even smaller, we may neglect all but the first terms in the above and so arrive at 
€, = 
To take a concrete example, if LZ, is one meter and x; and x2 are each one centimeter, then x2 is one hundred 


2 2 
times ty +02" 


With regard to the second spring, arguing as above, its elongation is (approximately) its stretch along 
its initial direction. As its initial direction is horizontal, its elongation is just the difference of the respective 
horizontal end displacements, namely, 


€2 = %3—- Wy 


Finally, the elongation of the third spring is (approximately) the difference of its respective vertical end 
displacements, i.e., 


€3 = 4 


We encode these three elongations in 


0 100 
e=Azx where A= -1 0410 
0 00 1 


Hooke’s law (Definition: "Hooke’s Law", p. 20) is an elemental piece of physics and is not perturbed by 
our leap from uniaxial to biaxial structures. The upshot is that the restoring force in each spring is still 
proportional to its elongation, i.e., y; = kje; where k; is the stiffness of the jth spring. In matrix terms, 
ky O O 
y=Ke where K= 0 k. O 
0 O kg 
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Balancing horizontal and vertical forces at mj, brings 


(—y2) — fi =0 
and 
yi — fo =0 
while balancing horizontal and vertical forces at mg brings 
y2— fs =0 
and 
y3 — fg =0 
We assemble these into 
0 -1 O 
1 0 0 
By=f where B= ; 
0 1 O 
0 O 1 


and recognize, as expected, that B is nothing more than A’. Putting the pieces together, we find that x 
must satisfy Sx = f where 


ee. OS ee 
0 k oO 0 
AAS ‘i 
a ko 


Applying one step of Gaussian Elimination (Section 3.1.3: Gaussian Elimination and the Uniaxial Truss) 
brings 


ko 0 —ko 0 “Ly fi 
0 ky 0 7) = fo 
0 0 v3 fit fs 
0 0 kg LA fa 
and back substitution delivers F 
a. J4 
w4 = kg 
0=fitfs 
ey aide 
2 i 
1-23 = fi 
x 3 Ks 


The second of these is remarkable in that it contains no components of x. Instead, it provides a condition 
on f. In mechanical terms, it states that there can be no equilibrium unless the horizontal forces on the two 
masses are equal and opposite. Of course one could have observed this directly from the layout of the truss. 
In modern, three—dimensional structures with thousands of members meant to shelter or convey humans 
one should not however be satisfied with the ‘visual’ integrity of the structure. In particular, one desires 
a detailed description of all loads that can, and, especially, all loads that can not, be equilibrated by the 
proposed truss. In algebraic terms, given a matrix S,, one desires a characterization of 
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1. all those f for which Sx = f possesses a solution 
2. all those f for which Sa = f does not possess a solution 


We will eventually provide such a characterization in our later discussion of the column space (Section 4.1) 
of a matrix. 

Supposing now that f; + f3 = 0 we note that although the system above is consistent it still fails to 
uniquely determine the four components of x. In particular, it specifies only the difference between x; and 
v3. As a result both 


p 0 
fa fo 
i Ra and «= “ 
0 a 
fa fa 
kg k3 


(3.10) 


XR 
Ill 
oF Oo Ff 


and still have a solution of Sx = f. Searching for the source of this lack of uniqueness we observe some 

redundancies in the columns of S. In particular, the third is simply the opposite of the first. As S is simply 
A’ K A we recognize that the original fault lies with A, where again, the first and third columns are opposites. 
These redundancies are encoded in z in the sense that 


Interpreting this in mechanical terms, we view z as a displacement and Az as the resulting elongation. In 
Az =0 


we see a nonzero displacement producing zero elongation. One says in this case that the truss deforms 
without doing any work and speaks of z as an unstable mode. Again, this mode could have been observed 
by a simple glance at Figure 3.2 (A simple swing). Such is not the case for more complex structures and so 
the engineer seeks a systematic means by which all unstable modes may be identified. We shall see later 
that all these modes are captured by the null space (Section 4.2) of A. 
From 
Sz=0 


one easily deduces that S$ is singular (Definition: "singular matrix", p. 23). More precisely, if S~' were to 
exist then S~1Sz would equal S~!0, i.e., z = 0, contrary to (3.10). As a result, Matlab will fail to solve 
Sa = f even when f is a force that the truss can equilibrate. One way out is to use the pseudo-inverse, 
as we shall see in the General Planar Truss (Section 3.3) module. 


3.3 The General Planar Truss’ 


Let us now consider something that resembles the mechanical prospection problem introduced in the intro- 
duction to matrix methods for mechanical systems (Section 3.1.1: Introduction). In the figure below we offer 
a crude mechanical model of a planar tissue, say, e.g., an excised sample of the wall of a vein. 


“This content is available online at <http://cnx.org/content/m10148/2.9/>. 
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A crude tissue model 


1 \ 
[> } g 20 ig 
10 
13 |14 18 
11 7 
4 L 15 
1 
. \5 } ” 6 
1 9 
2 4 . : 
1 1 6 L 
i 3 \ \ 
ea ied : .3 } 
Figure 3.3 


Elastic fibers, numbered 1 through 20, meet at nodes, numbered 1 through 9. We limit our observation to 
the motion of the nodes by denoting the horizontal and vertical displacements of node j by x2; (horizontal) 
and x2; (vertical), respectively. Retaining the convention that down and right are positive we note that the 
elongation of fiber 1 is 

€, = %— Xs 


while that of fiber 3 is 
€3 = %3 — 4. 


As fibers 2 and 4 are neither vertical nor horizontal their elongations, in terms of nodal displacements, are 
not so easy to read off. This is more a nuisance than an obstacle however, for noting our discussion of 
elongation (p. 25) in the small planar truss module, the elongation is approximately just the stretch along 
its undeformed axis. With respect to fiber 2, as it makes the angle — 7 with respect to the positive horizontal 
axis, we find 


7 ‘ as Lg — X1 + XQ — L190 
€2 = XgCOS ~ a — @19S1n = 


4 22 
Similarly, as fiber 4 makes the angle —2 with respect to the positive horizontal axis, its elongation is 


oe 3a . 37 3 —-X%7+%4— 2g 
4 = 2 —— } — xgsin = ; 
enrn 4 : A Noy) 


These are both direct applications of the general formula 
€j = L2n—1COS (0;) aa L2nSin (0;) (3.11) 


for fiber 7, as depicted in Figure 3.4 below, connecting node m to node n and making the angle 0; with 

the positive horizontal axis when node m is assumed to lie at the point (0,0). The reader should check that 
our expressions for e; and e3 indeed conform to this general formula and that eg and e4 agree with ones 
intuition. For example, visual inspection of the specimen suggests that fiber 2 can not be supposed to stretch 
(i.e., have positive e2) unless xg > 2 and/or x2 > £19. Does this jive with (3.11)? 
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original 


deformed 


Figure 3.4: Elongation of a generic bar, see (3.11). 


Applying (3.11) to each of the remaining fibers we arrive at e = Ax where A is 20-by-18, one row for 
each fiber, and one column for each degree of freedom. For systems of such size with such a well defined 
structure one naturally hopes to automate the construction. We have done just that in the accompanying 
M-file® and diary® . The M-file begins with a matrix of raw data that anyone with a protractor could have 
keyed in directly from Figure 3.3 (A crude tissue model): 


data = [ % one row of data for each fiber, the 

1 4 -pi/2 % first two columns are starting and ending 

1 5 -pi/4 % node numbers, respectively, while the third is the 
1 2 0 % angle the fiber makes with the positive horizontal 
2 4 -3*pi/4 

...and so on... a] 


This data is precisely what (3.11) requires in order to know which columns of A receive the proper cos or 
sin. The final A matrix is displayed in the diary” . 
The next two steps are now familiar. If kK denotes the diagonal matrix of fiber stiffnesses and f denotes 
the vector of nodal forces then 
y = Ke and ATy= f 


and so one must solve Sa = f where S = A’ KA. In this case there is an entire three-dimensional class of 
z for which Az = 0 and therefore Sz = 0. The three indicates that there are three independent unstable 
modes of the specimen, e.g., two translations and a rotation. As a result S is singular and x = S\f in 
MATLAB will get us nowhere. The way out is to recognize that S has 18 — 3 = 15 stable modes and that 
if we restrict S to ’act’ only in these directions then it ‘should’ be invertible. We will begin to make these 
notions precise in discussions on the Fundamental Theorem of Linear Algebra. (Section 4.5) For now let us 


®http://cnx.rice.edu/modules/m10148/latest/fiber.m 
Shttp://cnx.rice.edu/modules/m10148 /latest /lec2adj 
“http://cnx.rice.edu/modules/m10148 /latest /lec2adj 


axis 
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note that every matrix possesses such a pseudo-inverse and that it may be computed in MATLAB via the 
pinv command. Supposing the fiber stiffnesses to each be one and the edge traction to be of the form 


T 
f=(-110111-100010-1-10-11 ae 


we arrive at x via x=pinv(S)*f and offer below its graphical representation. 


3.3.1 Before-After Plot 


Figure 3.5: Before and after shots of the truss in Figure 3.3 (A crude tissue model). The solid (dashed) 
circles correspond to the nodal positions before (after) the application of the traction force, f. 


3.4 CAAM 335 Chapter 2 Exercises® 


Exercise 3.1 
With regard to the uniaxial truss figure (Figure 3.1: A uniaxial truss), 


e (i) Derive the A and K matrices resulting from the removal of the fourth spring, 

e (ii) Compute the inverse, by hand via Gauss-Jordan (Section 3.1.5: Gauss-Jordan Method: 
Computing the Inverse of a Matrix), of the resulting A’ KA with ky = ky = k3 =k 

e (iii) Use the result of (ii) to find the displacement. corresponding to the load f = (0,0, F)’. 


Exercise 3.2 

Generalize example 3, the general planar truss (Section 3.3), to the case of 16 nodes connected by 
42 fibers. Introduce one stiff (say k = 100) fiber and show how to detect it by ’properly’ choosing f. 
Submit your well-documented M-file as well as the plots, similar to the before-after plot (Figure 3.6) 
in the general planar module (Section 3.3), from which you conclude the presence of a stiff fiber. 


’This content is available online at <http://cnx.org/content /m10300/2.6/>. 


Figure 3.6: 


A copy of the before-after figure from the general planar module. 
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Chapter 4 


The Fundamental Subspaces 


4.1 Column Space’ 


4.1.1 The Column Space 


We begin with the simple geometric interpretation of matrix-vector multiplication. Namely, the multipli- 
cation of the n-by-1 vector x by the m-by-n matrix A produces a linear combination of the columns of A. 
More precisely, if a; denotes the jth column of A, then 


v2 
Ar = ( Qi 6-12. “Gz ) 


I 


©1Qa, + £202 +++: +LnAn 


The picture that I wish to place in your mind’s eye is that Az lies in the subspace spanned (Definition: 
"Span", p. 44) by the columns of A. This subspace occurs so frequently that we find it useful to distinguish 
it with a definition. 


Definition 4.1: Column Space 

The column space of the m-by-n matrix S is simply the span of the its columns, i.e. Ra(S) = 
{Sx |a¢€R”}. This is a subspace (Section 4.7.2) of R™. The notation Ra stands for range in this 
context. 


4.1.2 Example 


Let us examine the matrix: 


0 10 0 
A=|/ -1 01 0 (4.2) 
QAO) 0 4 


'This content is available online at <http://cnx.org/content /m10266/2.9/>. 
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The column space of this matrix is: 


0 1 0 0 
Ra(A)=¢ 21] -1 | +22] 0 | +23]| 1 | +aa| 0 | | ze R* (4.3) 
0 0 0 1 


Ra(A)=¢ 21] 1 | +22] 0 | +23] 0 | |xeR® (4.4) 


As the three remaining columns are linearly independent (Definition: "Linear Independence", p. 44) we 
may go no further. In this case, Ra (A) comprises all of R°. 


4.1.3 Method for Finding a Basis 


To determine the basis for Ra(A) (where A is an arbitrary matrix) we must find a way to discard its 
dependent columns. In the example above, it was easy to see that columns 1 and 3 were colinear. We seek, 
of course, a more systematic means of uncovering these, and perhaps other less obvious, dependencies. Such 
dependencies are more easily discerned from the row reduced form. In the reduction of the above problem, 
we come very easily to the matrix 


1 0 
0 0 (4.5) 
1 


Once we have done this, we can recognize that the pivot column are the linearly independent columns of 
Ayea- One now asks how this might help us distinguish the independent columns of A. For, although the 
rows of A;eg are linear combinations of the rows of A, no such thing is true with respect to the columns. 
The answer is: pay attention to the indices of the pivot columns. In our example, columns {1, 2, 4} 
are the pivot columns of A;eq and hence the first, second, and fourth columns of A, i.e., 


0 1 0 
—1 ’ 0 ’ 0 (4.6) 
0 0 1 


comprise a basis for Ra (A). In general: 


Definition 4.2: A Basis for the Column Space 
Suppose A is m-by-n. If columns {c; | j =1,...,7} are the pivot columns of A;eq then columns 
{c; | j =1,...,r} of A constitute a basis for Ra (A). 


4.2 Null Space’ 


4.2.1 Null Space 


Definition 4.3: Null Space 
The null space of an m-by-n matrix A is the collection of those vectors in R” that A maps to the 


?This content is available online at <http://cnx.org/content /m10293/2.9/>. 
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zero vector in R™. More precisely, 


NW (A) ={x eR" | Ax =0} 


4.2.2 Null Space Example 


As an example, we examine the matrix A: 


It is fairly easy to see that the null space of this matrix is: 


N(A)=ht |tER (4.8) 


oO FPF OO Fe 


This is a line in R*. 

The null space answers the question of uniqueness of solutions to Sx = f. For, if Sx = f and Sy = f 
then S (a —y) = Sa — Sy = f — f =O and so (4 — y) € -V (S). Hence, a solution to Sx = f will be unique 
if, and only if, WY (S) = {0}. 


4.2.3 Method for Finding the Basis 


Let us now exhibit a basis for the null space of an arbitrary matrix A. We note that to solve Ax = 0 is to 
solve A;eg® = 0. With respect to the latter, we suppose that 


{ej |F={1,...,r}} (4.9) 
are the indices of the pivot columns and that 
{ce |j={r+l,...,n}} (4.10) 
are the indices of the nonpivot columns. We accordingly define the r pivot variables 
eee a Cy oe (4.11) 
and the n —r free variables 


|e ees ee at Pewee 3 a (4.12) 


One solves A;eg2 = 0 by expressing each of the pivot variables in terms of the nonpivot, or free, variables. 
In the example above, x1, £2, and x4 are pivot while x3 is free. Solving for the pivot in terms of the free, we 
find x4 = 0, x3 = 21, x2 = 0, or, written as a vector, 


(4.13) 


8 
lI 
8 
w 
oOo FF OO Fe 
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where 73 is free. As 23 ranges over all real numbers the x above traces out a line in R+. This line is precisely 
the null space of A. Abstracting these calculations we arrive at: 


Definition 4.4: A Basis for the Null Space 

Suppose that A is m-by-n with pivot indices {c;|j={1,...,r}} and free indices 
{e; |g ={r+1,...,n}}. A basis for VW (A) may be constructed of n—r vectors {z',z?,...,2"-"} 
where z*, and only z*, possesses a nonzero in its c,4, component. 


4.2.4 A MATLAB Observation 


As usual, MATLAB has a way to make our lives simpler. If you have defined a matrix A and want to 
find a basis for its null space, simply call the function null (A). One small note about this function: if one 
adds an extra flag, ’r’, as in null(A, ’r’), then the basis is displayed "rationally" as opposed to purely 
mathematically. The MATLAB help pages define the difference between the two modes as the rational mode 
being useful pedagogically and the mathematical mode of more value (gasp!) mathematically. 


4.2.5 Final thoughts on null spaces 


There is a great deal more to finding null spaces; enough, in fact, to warrant another module. One important 
aspect and use of null spaces is their ability to inform us about the uniqueness of solutions. If we use the 
column space* to determine the existence of a solution x to the equation Ar = b. Once we know that a 
solution exists it is a perfectly reasonable question to want to know whether or not this solution is the only 
solution to this problem. The hard and fast rule is that a solution x is unique if and only if the null space of 
A is empty. One way to think about this is to consider that if Ax = 0 does not have a unique solution then, 
by linearity, neither does Ax = b. Conversely, if (Az =0) and (z#0) and (Ay=b) then A(z+y)=b 
as well. 


4.3 The Null and Column Spaces: An Example’ 


4.3.1 Preliminary Information 


Let us compute bases for the null and column spaces of the adjacency matrix associated with the ladder 
below. 


3 <http://cnx.org/content/columnspace/latest /> 
“This content is available online at <http://cnx.org/content /m10368/2.4/>. 
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An unstable ladder? 


Figure 4.1 


The ladder has 8 bars and 4 nodes, so 8 degrees of freedom. Denoting the horizontal and vertical 
displacements of node j by x2;-1 and x;, respectively, we arrive at the A matrix 


1 0 0 O 0 0 0 0 
-1 0 1 0 0 0 0 0 
0 O -1 0 0 0 0 0 
ge: 0 -1 0 0 0 1 0 0 
0 O O -!1 0 0 0 #1 
0 O O 0 1 0 0 0O 
0 O 0 -1 0 1 0O 
0 O 0 0 O -1 O 


4.3.2 Finding a Basis for the Column Space 


To determine a basis for R (A) we must find a way to discard its dependent columns. A moment’s reflection 
reveals that columns 2 and 6 are colinear, as are columns 4 and 8. We seek, of course, a more systematic 
means of uncovering these and perhaps other less obvious dependencies. Such dependencies are more easily 
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discerned from the row reduced form (Section 4.7.3) 


100 00 0 0 0 
0 10 0 0 -1 0 =O 
00100 0 0 +0 
Aue eer aye 0001 0 0 0 -1 
00001 0 0 +0 
00000 0 1 +0 
00000 0 0 =O 
00000 0 0 =O 


Recall that rref performs the elementary row operations necessary to eliminate all nonzeros below the 
diagonal. For those who can’t stand to miss any of the action I recommend rrefmovie. 

Each nonzero row of A;eg is called a pivot row. The first nonzero in each row of Ayeq is called a pivot. 
Each column that contains a pivot is called a pivot column. On account of the staircase nature of A,.q we 
find that there are as many pivot columns as there are pivot rows. In our example there are six of each and, 
again on account of the staircase nature, the pivot columns are the linearly independent (Definition: "Linear 
Independence", p. 44) columns of Ayeq. One now asks how this might help us distinguish the independent 
columns of A. For, although the rows of A;eq are linear combinations of the rows of A, no such thing is true 
with respect to the columns. In our example, columns {1,2,3,4,5, 7} are the pivot columns. In general: 


Proposition 4.1: 

Suppose A is m-by-n. If columns {c; | j =1,...,r} are the pivot columns of A;eq. then columns 
{cj |j =1,..,r} of A constitute a basis for R (A). 

Proof: 

Note that the pivot columns of A;eg. are, by construction, linearly independent. Suppose, however, 
that columns {c; | j =1,...,r} of A are linearly dependent. In this case there exists a nonzero 
xz € R” for which Ax = 0 and 


m=0, k¢{e;|f=1,.,r} (4.14) 


Now Az = 0 necessarily implies that A,.qz = 0, contrary to the fact that columns {c; | j = 1,...,r} 
are the pivot columns of Aye. 

We now show that the span of columns {c; | j =1,...,r} of A indeed coincides with R (A). 
This is obvious if r = n, i.e., if all of the columns are linearly independent. If r < n, there exists a 
q¢i{ce|j=1,...,r}. Looking back at Arq. we note that its gth column is a linear combination 
of the pivot columns with indices not exceeding g. Hence, there exists an x satisfying (4.14) and 
Areat = 0, and a, = 1. This x then necessarily satisfies Ax = 0. This states that the qth column 
of A is a linear combination of columns {c; | j =1,...,r} of A. 


4.3.3 Finding a Basis for the Null Space 


Let us now exhibit a basis for (A). We exploit the already mentioned fact that 1 (A) = -V (Area). 
Regarding the latter, we partition the elements of x into so called pivot variables, 


ee || oP acai 


and free variables 


{2% | k€ {cj | fj =1,...,r}} 
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There are evidently n — r free variables. For convenience, let us denote these in the future by 
eee | j =P bcm } 


One solves Ayeqx = 0, by expressing each of the pivot variables in terms of the nonpivot, or free, variables. 
In the example above, £1, £2, £3, £4, V5, and x7 are pivot while xg and 2g are free. Solving for the pivot in 
terms of the free we find 


r7 =0 
t5 = 0 
L4 = XB 
x3 = 0 
v2 = XE 
r1=0 


or, written as a vector, 


+ Xg (4.15) 


8 

II 

8 

fo) 
ScoOoOrFr coo OF oO 
FOoOoOoCocoHrFHOC Oo 


where xg and 2g are free. As xg and xg range over all real numbers, the x above traces out a plane in R®. 

This plane is precisely the null space of A and (4.15) describes a generic element as the linear combination 
of two basis vectors. Compare this to what MATLAB returns when faced with nul1(A,’r’). Abstracting 
these calculations we arrive at 


Proposition 4.2: 

Suppose that A is m-by-n with pivot indices {c;|j=1,...,r} and free indices 
{c;|j=r+1,...,n}. A basis for 4 (A). may be constructed of n —r vectors {z1, 2?,...,2”-"} 
where z*, and only z*, possesses a nonzero in its c,4, component. 


4.3.4 The Physical Meaning of Our Calculations 


Let us not end on an abstract note however. We ask what R(A) and -¥ (A) tell us about the ladder. 
Regarding R (A) the answer will come in the next chapter. The null space calculation however has revealed 
two independent motions against which the ladder does no work! Do you see that. the two vectors in (4.15) 
encode rigid vertical motions of bars 4 and 5 respectively? As each of these lies in the null space of A, the 
associated elongation is zero. Can you square this with the ladder as pictured in Figure 4.1 (An unstable 
ladder?)? I hope not, for vertical motion of bar 4 must ’stretch’ bars 1, 2, 6, and 7. How does one resolve 
this (apparent) contradiction? 
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4.4 Left Null Space’ 
4.4.1 Left Null Space 


If one understands the concept of a null space (Section 4.2), the left null space is extremely easy to understand. 


Definition 4.5: Left Null Space 
The Left Null Space of a matrix is the null space (Section 4.2) of its transpose, i.e., 


W (AT) = {yeER”" | Ay =0} 


The word "left" in this context stems from the fact that A?y = 0 is equivalent to y’A = 0 where y 
"acts" on A from the left. 
4.4.2 Example 


As Area was the key to identifying the null space (Section 4.2) of A, we shall see that A’, is the key to the 
null space of A’. If 


1 1 
AS | 4S (4.16) 
1 3 
then 
111 
Ab= (4.17) 
1 2 3 
and so 
11 
AS (4.18) 
0 1 2 
We solve AZ, = 0 by recognizing that y, and yg are pivot variables while y3 is free. Solving AZ,,y = 0 for 
the pivot in terms of the free we find yg = — (2y3) and y; = ys hence 
1 
N (AT) =4 ys| -2 | lyseR (4.19) 
1 


4.4.3 Finding a Basis for the Left Null Space 


The procedure is no different than that used to compute the null space (Section 4.2) of A itself. In fact 
Definition 4.6: A Basis for the Left Null Space 
Suppose that A’ is n-by-m with pivot indices {c,;|j={1,...,r}} and free indices 
{cj |j={r+l,...,m}}. A basis for 4 (A™) may be constructed of m — r_ vectors 
{z1,2?,...,2™-"} where z*, and only z*, possesses a nonzero in its cy4 Component. 


>This content is available online at <http://cnx.org/content /m10292/2.7/>. 
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4.5 Row Space* 
4.5.1 The Row Space 


As the columns of A? are simply the rows of A we call Ra (AT) the row space of A? . More precisely 


Definition 4.7: Row Space 
The row space of the m-by-n matrix A is simply the span of its rows, i.e., 


Ra(A™) = {ATy|yeR™} 


This is a subspace of R”. 


4.5.2 Example 


Let us examine the matrix: 


0 1 0 0 
A=] -10 10 (4.20) 
0 00 1 
The row space of this matrix is: 
0 -1 0 
1 0 0 
RA(AT) = 4 m1 + y2 + ys lye R? (4.21) 
0 1 0 
0 0 1 


As these three rows are linearly independent (Definition: "Linear Independence", p. 44) we may go no 
further. We "recognize" then RO (A™) as a three dimensional subspace (Section 4.7.2) of R*. 


4.5.3 Method for Finding the Basis of the Row Space 
Regarding a basis for RO (AT) we recall that the rows of A;ea, the row reduced form (Section 4.7.3) of the 
matrix A, are merely linear combinations of the rows of A and hence 


RA(A™) = RA (Area) (4.22) 


This leads immediately to: 


Definition 4.8: A Basis for the Row Space 
Suppose A is m-by-n. The pivot rows of Areq constitute a basis for RA (A’). 


With respect to our example, 


(4.23) 


Co Oo fF O&O 
e- © 
ee Oo Oo O&O 


comprises a basis for RA (AT). 
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4.6 Exercises: Columns and Null Spaces’ 


Exercises 


1. I encourage you to use rref and null for the following. 


e (i) Add a diagonal crossbar between nodes 3 and 2 in the unstable ladder figure (Figure 4.1: An 
unstable ladder?) and compute bases for the column and null spaces of the new adjacency matrix. 
As this crossbar fails to stabilize the ladder, we shall add one more bar. 
e (ii) To the 9 bar ladder of (i) add a diagonal cross bar between nodes 1 and the left end of bar 6. 
Compute bases for the column and null spaces of the new adjacency matrix. 
2. We wish to show that N (A) = N (A™ A) regardless of A. 
e (i) We first take a concrete example. Report the findings of null when applied to A and A’ A 
for the A matrix associated with the unstable ladder figure (Figure 4.1: An unstable ladder?). 
e (ii) Show that N (A) C N (A7A), ie. that if Az = 0 then A? Az = 0. 
e (iii) Show that N (A7A) C N (A), ie., that if ATA =0 then Ax = 0. (Hint: if ATAzx = 0 then 
x? AT Ar = 0.) 
3. Suppose that A is m-by-n and that N (A) = R”. Argue that A must be the zero matrix. 


4.7 Appendices /Supplements 


4.7.1 Vector Space® 
4.7.1.1 Introduction 


You have long taken for granted the fact that the set of real numbers, R, is closed under addition and 
multiplication, that each number has a unique additive inverse, and that the commutative, associative, and 
distributive laws were right as rain. The set, C, of complex numbers also enjoys each of these properties, as 
do the sets R” and C” of columns of n real and complex numbers, respectively. 

To be more precise, we write xz and y in R” as 

uv = (#1, %2,..., Ln) 


T 
Y = (Y1,Y2,-++5 Yn) 
and define their vector sum as the elementwise sum 


TTY 
T+ Y2 

ety = ’ (4.24) 
In Yn 


and similarly, the product of a complex scalar, z € C with 2 as: 


20 = : (4.25) 


“This content is available online at <http://cnx.org/content /m10367/2.4/>. 
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4.7.1.2 Vector Space 
These notions lead naturally to the concept of vector space. A set V is said to be a vector space if 


x+y=y-+-2 for each x and y in V 

at+y+z2=a2+y+4 <2 for each z, y, and z in V 

There is a unique "zero vector" such that «+0 =< for each x in V 
For each x in V there is a unique vector —x such that 2 + —x = 0. 
la=2 

(c1C2) & = cy (cgx) for each x in V and c; and cg in C. 

c(a+y) = cx + cy for each x and y in V and c in C. 

(cy +2) u = cx + Cox for each x in V and c, and cg in C. 


OOS Oe Ne 


4.7.2 Subspaces’ 
4.7.2.1 Subspace 


A subspace is a subset of a vector space (Section 4.7.1) that is itself a vector space. The simplest example 
is a line through the origin in the plane. For the line is definitely a subset and if we add any two vectors on 
the line we remain on the line and if we multiply any vector on the line by a scalar we remain on the line. 
The same could be said for a line or plane through the origin in 3 space. As we shall be travelling in spaces 
with many many dimensions it pays to have a general definition. 


Definition 4.9: A subset S of a vector space V is a subspace of V when 
1. if 2 and y belong to S then so does x + y. 
2. if x belongs to S and t is real then tx belong to S. 
As these are oftentimes unwieldy objects it pays to look for a handful of vectors from which the entire 
subset may be generated. For example, the set of x for which x; + 72 + ©3 + x4 = 0 constitutes a subspace 
of R*. Can you ’see’ this set? Do you ’see’ that 


and 


and 


not only belong to a set but in fact generate all possible elements? More precisely, we say that these vectors 
span the subspace of all possible solutions. 


®This content is available online at <http://cnx.org/content /m10297/2.6/>. 
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Definition 4.10: Span 
A finite collection {s1, 82,...,,} of vectors in the subspace S' is said to span S' if each element of 
S can be written as a linear combination of these vectors. That is, if for each s € S there exist n 
reals {x1,%2,...,2%n} such that s = 7181 + g82 +--+ + an Sn. 
When attempting to generate a subspace as the span of a handful of vectors it is natural to ask what is 
the fewest number possible. The notion of linear independence helps us clarify this issue. 


Definition 4.11: Linear Independence 

A finite collection {s1, 52,...,8,} of vectors is said to be linearly independent when the only reals, 
{21,@2,...,@,} for which 71 +22 +---+ a, = 0 are v1 = 4g =--+ = &, = 0. In other words, when 
the null space (Section 4.2) of the matrix whose columns are {81, 52,...,5,} contains only the zero 
vector. 


Combining these definitions, we arrive at the precise notion of a ’generating set.’ 


Definition 4.12: Basis 
Any linearly independent spanning set of a subspace S is called a basis of S. 


Though a subspace may have many bases they all have one thing in common: 


Definition 4.13: Dimension 
The dimension of a subspace is the number of elements in its basis. 


4.7.3 Row Reduced Form” 
4.7.3.1 Row Reduction 


A central goal of science and engineering is to reduce the complexity of a model without sacrificing its 
integrity. Applied to matrices, this goal suggests that we attempt to eliminate nonzero elements and so 
*uncouple’ the rows. In order to retain its integrity the elimination must obey two simple rules. 

Definition 4.14: Elementary Row Operations 

1. You may swap any two rows. 

2. You may add to a row a constant multiple of another row. 

With these two elementary operations one can systematically eliminate all nonzeros below the diagonal. 

For example, given 


0 1 0 0 
-1 01 0 
(4.26) 
0 0 0 1 
1 23 4 
it seems wise to swap the first and fourth rows and so arrive at 
1 23 4 
0 1 0 0 
(4.27) 
-1 01 0 
0 0 0 1 
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adding the first row to the third now produces 


123 4 
0 10 0 
(4.28) 
02 44 
0 0 0 1 
subtracting twice the second row from the third yields 
(4.29) 


3.4 
0 0 
4 4 
0 1 


eo Oo Oo F 
Co oO fF W& 


a matrix with zeros below its diagonal. This procedure is not restricted to square matrices. For example, 
given 


1 
2 
3 


o eR 
o -. Re 


i 
2 (4.30) 
3 


we start at the bottom left then move up and right. Namely, we subtract 3 times the first row from the 
third and arrive at 


1111 
24 4 2 (4.31) 
0 2 2 0 

and then subtract twice the first row from the second, 
TS hy de pad 
0 2 2 0 (4.32) 
0 2 2 0 

and finally subtract the second row from the third, 
11411 
0 2 2 0 (4.33) 
0 0 0 0 


It helps to label the before and after matrices. 
Definition 4.15: The Row Reduced Form 
Given the matrix A we apply elementary row operations until each nonzero below the diagonal is 
eliminated. We refer to the resulting matrix as Ayea. 
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4.7.3.2 Uniqueness and Pivots 


As there is a certain amount of flexibility in how one carries out the reduction it must be admitted that 
the reduced form is not unique. That is, two people may begin with the same matrix yet arrive at different 
reduced forms. The differences however are minor, for both will have the same number of nonzero rows and 
the nonzeros along the diagonal will follow the same pattern. We capture this pattern with the following 
suite of definitions, 

Definition 4.16: Pivot Row 

Each nonzero row of Areq is called a pivot row. 

Definition 4.17: Pivot 

The first nonzero term in each row of A;eq is called a pivot. 

Definition 4.18: Pivot Column 

Each column of A;eq that contains a pivot is called a pivot column. 

Definition 4.19: Rank 

The number of pivots in a matrix is called the rank of that matrix. 


Regarding our example matrices, the first (4.26) has rank 4 and the second (4.30) has rank 2. 


4.7.3.3 Row Reduction in MATLAB 


MATLAB’s rref command goes full-tilt and attempts to eliminate ALL off diagonal terms and to leave 
nothing but ones on the diagonal. I recommend you try it on our two examples. You can watch its individual 
decisions by using rrefmovie instead. 


Chapter 5 


Least Squares 


5.1 Least Squares’ 


5.1.1 Introduction 


We learned in the previous chapter that Ax = b need not possess a solution when the number of rows of A 
exceeds its rank, i.e., r << m. As this situation arises quite often in practice, typically in the guise of ’more 
equations than unknowns,’ we establish a rationale for the absurdity Ax = b. 


5.1.2 The Normal Equations 
The goal is to choose x such that Az is as close as possible to b. Measuring closeness in terms of the sum of 
the squares of the components we arrive at the ‘least squares’ problem of minimizing 
res 

(|| Ax — b ||)? = (Ax —b)? (Aa — b) (5.1) 
over all x € R. The path to the solution is illuminated by the Fundamental Theorem. More precisely, we 
write b=br+bn , be €R(A) and by €N(A7) . On noting that (i) (At—bg) € R(A) , ce R” 
and (ii)IR (A) L N(A7) we arrive at the Pythagorean Theorem. 
Pythagorean Theorem 


norm? (Ax—b) = (|| Ax —be+bn |l)” 


: : (5.2) 
(|| Av — bp II)" + (Il by II) 


It is now clear from the Pythagorean Theorem (5.2: Pythagorean Theorem) that the best x is the one that 


satisfies 


As bry € R(A) this equation indeed possesses a solution. We have yet however to specify how one computes 
br given b. Although an explicit expression for br, the so called orthogonal projection of b onto R (A), 
in terms of A and 6 is within our grasp we shall, strictly speaking, not require it. To see this, let us note 
that if x satisfies the above equation (5.3) then 


Az—b = Ax—brt+bn 65.4) 


= —bn 
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As by is no more easily computed than bg you may claim that we are just going in circles. The practical’ 
information in the above equation (5.4) however is that (Ax — b) € A”, i.e, AT (Ax — b) = 0, ice., 


A? Ax = A™b (5.5) 


As ATDER (AT) regardless of b this system, often referred to as the normal equations, indeed has 
a solution. This solution is unique so long as the columns of A” A are linearly independent, i.e., so long 
as N(A7A) = {0}. Recalling Chapter 2, Exercise 2, we note that this is equivalent to N(A) = {0}. We 
summarize our findings in 

Theorem 5.1: 
The set of « € bg for which the misfit (|| Aa —b ||)” is smallest is composed of those x for which 
A’ Ax = A™b. There is always at least one such x. There is exactly one such « if N(A) = {0}. 


1 1 1 
As a concrete example, suppose with reference to Figure 5.1 that A= | 0 1 | andb=| 1 
0 0 1 


Figure 5.1: The decomposition of b. 
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As b # R(A) there is no x such that Ax = b. Indeed, (|| Ax — 6 ||)? = (v1 + v2 + —1)?+(x2 — 1)? +1 > 1, 


0 
with the minimum uniquely attained at 7 = , in agreement with the unique solution of the above 
1 
. T 1 1 r 1 ; ae 
equation (5.5), for A* A = and A*b = . We now recognize, a posteriori, that by = Ax = 
1 2 2 
1 
1 is the orthogonal projection of b onto the column space of A. 
0 


5.1.3 Applying Least Squares to the Biaxial Test Problem 


We shall formulate the identification of the 20 fiber stiffnesses in this previous figure (Figure 3.3: A crude 
tissue model), as a least squares problem. We envision loading, f, the 9 nodes and measuring the associated 
18 displacements, x. From knowledge of « and f we wish to infer the components of K = diag (k) where k 
is the vector of unknown fiber stiffnesses. The first step is to recognize that 


A’ KAz = f (5.6) 


may be written as 
Bk=f , B=A'diag (Az) (5.7) 


Though conceptually simple this is not of great use in practice, for B is 18-by-20 and hence the above 
equation (5.7) possesses many solutions. The way out is to compute k as the result of more than one 
experiment. We shall see that, for our small sample, 2 experiments will suffice. To be precise, we suppose 
that. x! is the displacement produced by loading f! while x? is the displacement produced by loading f?. 
A’ diag (Act) fz ; ; 

and f = This B is 36-by-20 
Al diag (Az?) Cie 
and so the system Bk = f is overdetermined and hence ripe for least squares. 
We proceed then to assemble B and f. We suppose f! and f? to correspond to horizontal and vertical 
stretching 


We then piggyback the associated pieces in B = 


T 
fs(-lo0o0010-100010-100010) (5.8) 


de 
PAO 10 1G 10 O00 0 <7 6-0 et) (5.9) 


respectively. For the purpose of our example we suppose that each k; = 1 except kg = 5. We assemble 
A’ K Aas in Chapter 2 and solve 
A’ K Ax) = fi (5.10) 


with the help of the pseudoinverse. In order to impart some ‘reality’ to this problem we taint each x? with 
10 percent noise prior to constructing B. Regarding 


B’ Bk = B' f (5.11) 


we note that Matlab solves this system when presented with k=B\f when B is rectangular. We have plotted 
the results of this procedure in the Figure 5.2. The stiff fiber is readily identified. 
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fiber stiffness 
to fh 


ht 


5 10 15 20 25 
fiber number 


Figure 5.2: Results of a successful biaxial test. 


5.1.4 Projections 

From an algebraic point of view (5.5))is an elegant reformulation of the least squares problem. Though easy 
to remember it unfortunately obscures the geometric content, suggested by the word ’projection,’ of (5.4). 
As projections arise frequently in many applications we pause here to develop them more carefully. With 


respect to the normal equations we note that if N(A) = {0} then 
x= (ATA) ATD (5.12) 
and so the orthogonal projection of 6 onto R (A) is: 
be = Ax 
. (5.13) 


Defining 
AP (5.14) 
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(5.13) takes the form bg = Pb. Commensurate with our notion of what a ’projection’ should be we 
expect that P map vectors not in R(A) onto R(A) while leaving vectors already in R(A) unscathed. More 
succinctly, we expect that Pbr = br, i.e., Pbr = Pbr. As the latter should hold for all b € R™ we expect 
that 

Pe=P (5.15) 


With respect to (5.14) we find that indeed 


1 


Pp? = A(ATA)*ATA(ATA) AT 
= A(ATA)~*AT (5.16) 
= P 


We also note that the P in (5.14) is symmetric. We dignify these properties through 


Definition 5.1: orthogonal projection 
A matrix P that satisfies P? = P is called a projection. A symmetric projection is called an 
orthogonal projection. 
We have taken some pains to motivate the use of the word ’projection.’ You may be wondering however 
what symmetry has to do with orthogonality. We explain this in terms of the tautology 


b= Pb—Ib (5.17) 


Now, if P is a projection then so too is J — P. Moreover, if P is symmetric then the dot product of b’s two 
constituents is 


(Pb)’ (I—P)b = bTPT(I—P)b 
a ei) (5.18) 
= b70b 
= 0 


ie, Pb is orthogonal to (I—P)b. As examples of a nonorthogonal projections we offer P = 
1 0 O 
=> 0 0 and I — P. Finally, let us note that the central formula, P = A(ATA)~* = A’, is 


even a bit more general than advertised. It has been billed as the orthogonal projection onto the column 
space of A. The need often arises however for the orthogonal projection onto some arbitrary subspace M. 
The key to using the old P is simply to realize that every subspace is the column space of some matrix. 
More precisely, if 

{£1, ..., Lm} (5.19) 


is a basis for M then clearly if these x; are placed into the columns of a matrix called A then R(A) = M. 


T 
For example, if M is the line through ( 1 1 ) then 


a 4(1 1) 
(5.20) 
=. aie 4 
2, es ame 


is orthogonal projection onto M. 
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5.1.5 Exercises 


1. 


Gilbert Strang was stretched on a rack to lengths @ = 6, 7, and 8 feet under applied forces of f = 1, 2, 
and 4 tons. Assuming Hooke’s law @— L = cf, find his compliance, c, and original height, L, by least 
squares. 


. With regard to the example of §3 note that, due to the the random generation of the noise that taints 


the displacements, one gets a different ’answer’ every time the code is invoked. 


a. Write a loop that invokes the code a statistically significant number of times and submit bar plots 
of the average fiber stiffness and its standard deviation for each fiber, along with the associated 
M-file. 

b. Experiment with various noise levels with the goal of determining the level above which it becomes 
difficult to discern the stiff fiber. Carefully explain your findings. 


T 
. Find the matrix that projects R® onto the line spanned by ( 101 ) ; 


ze E. 
4. Find the matrix that projects R° onto the plane spanned by ( 101 ) and ( 11 -1 ) 


5. If P is the projection of R™ onto a k-dimensional subspace M, what is the rank of P and what is 


R(P)? 


Chapter 6 


Matrix Methods for Dynamical Systems 


6.1 Nerve Fibers and the Dynamic Strang Quartet’ 


6.1.1 Introduction 


Up to this point we have largely been concerned with 


1. Deriving linear systems of algebraic equations (from considerations of static equilibrium) and 
2. The solution of such systems via Gaussian elimination. 


In this module we hope to begin to persuade the reader that our tools extend in a natural fashion to the 
class of dynamic processes. More precisely, we shall argue that 


1. Matrix Algebra plays a central role in the derivation of mathematical models of dynamical systems 
and that, 

2. With the aid of the Laplace transform in an analytical setting or the Backward Euler method in the 
numerical setting, Gaussian elimination indeed produces the solution. 


6.1.2 Nerve Fibers and the Dynamic Strang Quartet 
6.1.2.1 Gathering Information 


A nerve fiber’s natural electrical stimulus is not direct current but rather a short burst of current, the so- 
called nervous impulse. In such a dynamic environment the cell’s membrane behaves not only like a leaky 
conductor but also like a charge separator, or capacitor. 


'This content is available online at <http://cnx.org/content /m10168/2.6/>. 
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An RC model of a nerve fiber 


Rj 
x; <. 
¥, 
C 2) fe C 
ydee JS" CP yt | 
Em %s 
Figure 6.1 


The typical value of a cell’s membrane capacitance is 


pF 
— 1 —_ 
. cm? 


where yF denotes micro-Farad. Recalling our variable conventions (Section 2.1.1: 


Strang Quartet), the capacitance of a single compartment is 


l 
Cy = rare 


Nerve Fibers and the 


and runs parallel to each R,,, see Figure 6.1 (An RC model of a nerve fiber). This figure also differs from 
the simpler circuit (Figure 2.3: The fully dressed circuit model) from the introductory electrical modeling 
module in that it possesses two edges to the left of the stimuli. These edges serve to mimic that portion of 
the stimulus current that is shunted by the cell body. If A, denotes the surface area of the cell body, then 


it has 


Definition 6.1: capacitance of cell body 
Cob = Ache 

Definition 6.2: resistance of cell body 
Reb = AcbPm- 


6.1.2.2 Updating the Strang Quartet 


We ask now how the static Strang Quartet (p. 8) of the introductory electrical module should be augmented. 


6.1.2.2.1 Updating (S1’) 


Regarding (S1’) (p. 13) we proceed as before. The voltage drops are 


€y = 


€6 = %2 — %3 


€7 = 73 


€g = 173 — Em 


and so 
0 -1 0 0O 
1 -1 0 0 
0 -1 1 0 
e=b—Ax where b=(-E,,) : and A= Se. oe 
1 0 -1 O 
0 0 -1 1 
0 0 oO -1 
1 0 oO -1 


6.1.2.2.2 Updating (S2) 
To update (S2) (Section 2.1.2.2: Strang Quartet, Step 2) we must now augment Ohm’s law with 


Definition 6.3: Voltage-current law obeyed by a capacitor 
The current through a capacitor is proportional to the time rate of change of the potential across 


it. 
This yields, (denoting derivative by ’), 


Y= Coper’ 
€2 
Y2>= 
Reb 
€3 
¥3 = R; 
te = On er 
€5 
Y5 = Rm 
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Y= Cme7’ 
€8 
Ys = Re, 
or, in matrix terms, 
y =Ge+Ce' 
where 
0 O 0 0 O 0 O O 
1 
Oo ae OF 30) 8 0 O O 
Ors 205, ee 80 OE 50 HO 0 
0 O 0 0 O 0 0 O 
G= je 
0 O OF Oh ge OP 0" 0 
0 O 0 0 O Be 0 0 
0 O 0 0 O 0 O O 
Oh 0. 0D 20 07 TO 
and 
Cp 0 0 0 0 0 0 0 
0 00 0 00 0 0 
0 00 0 00 0 0 
Cx 0 00 Cr 0 0 0 0 
0 00 0 00 0 0 
0 00 0 00 0 0 
0 00 0 00 GC, 0 
0 00 0 00 0 0 


are the conductance and capacitance matrices. 


6.1.2.2.3 Updating (S3) 


As Kirchhoff’s Current law is insensitive to the type of device occupying an edge, step (S3) proceeds exactly 
as before (Section 2.1.2.3: Strang Quartet, Step 3). 


to — Yi — Yo — y3 =O 


Y3 — Ya — Ys — Yo = O 


Yo — ¥7 — ys = 0 
or, in matrix terms, 
t0 
A’y=-—f where f=| 0 


57 


6.1.2.2.4 Step (S4): Assembling 
Step (S4) remains one of assembling, 

(ATy = —f) > (A (Ge + Ce’) = —f) = (AT (G(b- Az) + C (b’ — Az’)) = —f) 
becomes 


A’CAz' + ATGAx = A'Gb+ f + ATCH. (6.1) 


This is the general form of the potential equations for an RC circuit. It presumes of the user knowledge 
of the initial value of each of the potentials, 


x(0)=X (6.2) 
Regarding the circuit of Figure 6.1 (An RC model of a nerve fiber), and letting G = 7 we find 
Cy 0 0 Gop + Gi —G; 0 
ATCA=| 0 C 0], A’GA= -—G;i 2Gi+Gn = -G; 
0 0 C 0 —G; Git+Gm 
Geb 
APGd = Eye | Gy, and A’Cbv’=1| 0 
Gm 0 
and an initial (rest) potential of 
; T 
x(0)=Em} 1 
1 


6.1.2.3 Modes of Attack 


We shall now outline two modes of attack on such problems. The Laplace Transform (Section 6.2) is an 
analytical tool that produces exact, closed-form solutions for small tractable systems and therefore offers 
insight into how larger systems ’should’ behave. The Backward-Euler method (Section 6.4) is a technique 
for solving a discretized (and therefore approximate) version of (6.1). It is highly flexible, easy to code, and 
works on problems of great size. Both the Backward-Euler and Laplace Transform methods require, at their 
core, the algebraic solution of a linear system of equations. In deriving these methods we shall find it more 
convenient to proceed from the generic system 


ve =But+g (6.3) 
With respect to our fiber problem 
B = (-(A?CA)"") ATGA 
—(Gen+Gi) Gi 0 
Ce. Ce 
a G, ~(26:4G mn) Gi (6.4) 
m Cm Cus 
Gi (Gi Gm) 
0 Ce om 
and 
Ger Em tio 
= Cob 
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6.2 The Laplace Transform’ 


The Laplace Transform is typically credited with taking dynamical problems into static problems. Recall 
that the Laplace Transform of the function h is 


MATLAB is very adept at such things. For example: 
Example 6.1: The Laplace Transform in MATLAB 


> syms t 

> laplace(exp(t)) 
ans = 1/(s-1) 

> laplace(t*(exp(-t)) 
ans = 1/(st1)72 


The Laplace Transform of a matrix of functions is simply the matrix of Laplace transforms of the individual 
elements. 


Example 6.2: Laplace Transform of a matrix of functions 


t 1 
=I 
L = ‘ 
1 


—t 
te (+2 


Now, in preparing to apply the Laplace transform to our equation from the dynamic strang quartet module 
(6.3): 


x =Bxrt+g, 
we write it as 
dx 
£ ooo L(Ba+q) (6.5) 


and so must determine how £ acts on derivatives and sums. With respect to the latter it follows directly 
from the definition that 


L(Bx+q) L£(Ba)+ L(g) 


= BL(«)+Ll(qg) 


I 


(6.6) 


Regarding its effect. on the derivative we find, on integrating by parts, that 


dx ety Ax (t) o = 
—|\= (st) dt = x(t) e 90 + / (st). (t) dt. 
£(F) | € 7 x (t)e la +s € x (t) 


0 


Supposing that 2 and s are such that a (t)e~* — 0 as t > 00 we arrive at 


L (=) abs) 2 (0)) (6.7) 


?This content is available online at <http://cnx.org/content/m10169/2.5/>. 
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Now, upon substituting (6.6) and (6.7) into (6.5) we find 
s£(x) —2(0) = BL(2) + £(9), 
which is easily recognized to be a linear system for L(x), namely 
(sl — B) L(x) =L(g)+2(0). (6.8) 


The only thing that distinguishes this system from those encountered since our first brush (Section 2.1) with 
these systems is the presence of the complex variable s. This complicates the mechanical steps of Gaussian 
Elimination or the Gauss-Jordan Method but the methods indeed apply without change. Taking up the 
latter method, we write 
-1 
L(x) = (sI—B) (L(g) + 2(0)). 

The matrix (sI — B)"* is typically called the transfer function or resolvent, associated with B, at s. 
We turn to MATLAB for its symbolic calculation. (for more information, see the tutorial? on MATLAB’s 
symbolic toolbox). For example, 


Example 6.3 


> B= [2 -1; -1 2] 


V 
=e) 
I 


= inv(s*eye(2) -B) 


[ (s-2)/(s*s-4*s+3), -1/(s*s-4*s+3) ] 


[ -1/(s*s-4*s+3), (s-2)/(s*s-4*s+3) ] 


We note that (sI — B)* is well defined except at the roots of the quadratic, s*? — 4s + 3. This quadratic is 
the determinant of s/ — B and is often referred to as the characteristic polynomial of B. Its roots are 
called the eigenvalues of B. 


Example 6.4 
As a second example let us take the B matrix of the dynamic Strang quartet module (6.4) with 
the parameter choices specified in fib3.m* , namely 


-0.135 0.125 0 
B= 0.5 1.01 0.5 (6.9) 
0 0.5 —0.51 


The associated (sI — B)" is a bit bulky (please run fib3.m® ) so we display here only the denom- 
inator of each term, i.e., 
s? + 1.6558” + 0.40785 + 0.0039. (6.10) 


eee 
_ the 6 


Assuming a current stimulus of the form ig (t) = Fogg and Em = 0 brings 


0.191 
L(g) (s) = 0 


0 


3http://www.mathworks.com/access/helpdesk/help/toolbox/symbolic/symbolic.shtml 
“http://www.caam.rice.edu/~+caam335/cox/lectures/fib3.m 
Shttp://www.caam.rice.edu/~caam335/cox/lectures/fib3.m 
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and so (6.10) persists in 
s? + 1.58 + 0.27 
0.191 


£ (2) = (sI - BY L(g) = 7 0.58 + 0.26 
(s+) (s3 + 1.6555? + 0.4078s + 0.0039) ah 


Now comes the rub. A simple linear solve (or inversion) has left us with the Laplace transform of x. The 
accursed 


Theorem 6.1: No Free Lunch Theorem 
We shall have to do some work in order to recover x from L (2). 


confronts us. We shall face it down in the Inverse Laplace module (Section 6.3). 


6.3 The Inverse Laplace Transform’ 


6.3.1 To Come 
In The Transfer Function (Section 9.2) we shall establish that the inverse Laplace transform of a function h 
is 

1 


~ On 


L~' (h) (t) i. et¥Dth (c+ yi) t) dy (6.11) 


where i = \/2 — 1 and the real number c is chosen so that all of the singularities of h lie to the left of the 
line of integration. 


6.3.2 Proceeding with the Inverse Laplace Transform 


With the inverse Laplace transform one may express the solution of 2’ = Bxa+q , as 
x(t) = £7" ((s1 — B)"') (L(g) +2 (0)) (6.12) 
As an example, let us take the first component of £ (a), namely 


0.19 (s? + 1.5s + 0.27) 
(s + 4)* (s3 + 1.65552 + 0.40785 + 0.0039) 


Ley (s) = 


We define: 


Definition 6.4: poles 
Also called singularities, these are the points s at which L,, (s) blows up. 


These are clearly the roots of its denominator, namely 


273 
—1/100, — 329/400+ “-=, and —1/6. (6.13) 


All four being negative, it suffices to take c = 0 and so the integration in (6.11) proceeds up the imaginary 
axis. We don’t suppose the reader to have already encountered integration in the complex plane but hope 
that this example might provide the motivation necessary for a brief overview of such. Before that however 
we note that MATLAB has digested the calculus we wish to develop. Referring again to fib3.m’ for details 
we note that the ilaplace command produces 


°This content is available online at <http://cnx.org/content /m10170/2.8/>. 
“http://www.caam.rice.edu/~caam335/cox/lectures/fib3.m 
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(0.0554t3 + 4.5464t? + 1.085 +474.19)ee + 


) 


a (t) = 211.35e%m0  — 


“tw” (262.842cosh (FZ) ) + 262.836sinh ( 


140 


a J 


x (mV) 


60F J 


40 


IN 
pz 4 4 4 
i] §0 100 180 200 2m 300 
t (ms) 


Figure 6.2: The 3 potentials associated with the RC circuit model figure (Figure 6.1: An RC model 


of a nerve fiber). 


The other potentials, see the figure above, possess similar expressions. Please note that each of the poles 
of £(a,) appear as exponents in 2, and that the coefficients of the exponentials are polynomials whose 


degrees is determined by the order of the respective pole. 


6.4 The Backward-Euler Method* 


Where in the Inverse Laplace Transform (Section 6.3) module we tackled the derivative in 
xv =Brt+g, (6.14) 
via an integral transform we pursue in this section a much simpler strategy, namely, replace the derivative 
with a finite difference quotient. That is, one chooses a small dt and ’replaces’ (6.14) with 
z(t) -— x(t — dt z 
dt 
The utility of (6.15) is that it gives a means of solving for % at the present time, ¢, from the knowledge of 
& in the immediate past, t — dt. For example, as £ (0) = x (0) is supposed known we write (6.15) as 


(3-8) e(at) = 2 + 9 (ae) 


(6.15) 


dt 


’This content is available online at <http://cnx.org/content /m10171/2.6/>. 
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Solving this for % (dt) we return to (6.15) and find 


(5 = B) & (2dt) = ne) +9 (2dt). 


and solve for #(2dt). The general step from past to present, 


re (= 7 B) ee 16 (dt) (6.16) 


is repeated until some desired final time, Tdt, is reached. This equation has been implemented in fib3.m® 
with dt = 1 and B and g as in the dynamic Strang module (6.4). The resulting % ( run fib3.m!° yourself!) 
is indistinguishable from the plot we obtained (Figure 6.2) in the Inverse Laplace module. 

Comparing the two representations, this equation (6.12) and (6.16), we see that they both produce the 
solution to the general linear system of ordinary equations, see this eqn (6.3), by simply inverting a shifted 
copy of B. The former representation is hard but exact while the latter is easy but approximate. Of course 
we should expect the approximate solution, z , to approach the exact solution, x, as the time step, dt , 
approaches zero. To see this let us return to (6.16) and assume, for now, that g = 0 . In this case, one can 
reverse the above steps and arrive at the representation 


& (jdt) = (Z 2 apy)’ « (0) (6.17) 


Now, for a fixed time t we suppose that dt = ; and ask whether 


This limit, at least when B is one-by-one, yields the exponential 
x (t) = e?*x (0) 


clearly the correct solution to this equation (6.3). A careful explication of the matrix exponential and its 
relationship to this equation (6.12) will have to wait until we have mastered the inverse laplace transform. 


6.5 Exercises: Matrix Methods for Dynamical Systems” 


1. Compute, without the aid of a machine, the Laplace transforms of e' and te~’. Show ALL of your 

work. 

Extract from fib3.m analytical expressions for x2 and 73. 

3. Use eig to compute the eigenvalues of B as given in this equation (6.9). Use det to compute the 
characteristic polynomial of B. Use roots to compute the roots of this characteristic polynomial. 
Compare these to the results of eig. How does Matlab compute the roots of a polynomial? (type help 
roots for the answer). 

4. Adapt the Backward Euler portion of fib3.m so that one may specify an arbitrary number of com- 
partments, as in fib1.m. Submit your well documented M-file along with a plot of x; and x19 versus 
time (on the same well labeled graph) for a nine compartment fiber of length | = 1cm. 


ie) 


°http://www.caam.rice.edu/~caam335/cox/lectures/fib3.m 
Mhttp://www.caam.rice.edu/~caam335/cox/lectures/fib3.m 
11 This content is available online at <http://cnx.org/content /m10526/2.4/>. 
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5. Derive this equation (6.17) from a previous equation (6.16) by working backwards toward x (0). Along 
the way you should explain why 


: a5 
6. Show, for scalar B, that ((1 - +B) ) — ePt as j — oo. Hint: By definition 


t a ts jlog—— 
((: = “B) =e Mage 
J 


now use L’Hopital’s rule to show that jlog—s — Bt. 
j 


6.6 Supplemental 


6.6.1 Matrix Analysis of the Branched Dendrite Nerve Fiber™ 
6.6.1.1 Introduction 


In the prior modules on static (Section 2.1) and dynamic electrical systems, we analyzed basic, hypothetical 
one-branch nerve fibers using a modeling methodology we dubbed the Strang Quartet. You may be asking 
yourself whether this method is stout enough to handle the real fiber of our minds. Indeed, can we use our 
tools in a real-world setting? 


6.6.1.2 Presentation 


An Actual Nerve Fiber 


Figure 6.3: A pyramidal neuron from the CA3 region of a rat’s hippocampus, scanned at (FIX ME) X 
magnification. 


!2This content is available online at <http://cnx.org/content/m10177/2.7/>. 
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To answer your question, the above is a rendering of a neuron from a rat’s hippocampus. The tools we 
have refined will enable us to model the electrical properties of a dendrite leaving the neuron’s cell body. A 
three-branch model of such a dendrite, traced out with painstaking accuracy, appears in the diagram below. 


3-branch Dendrite Model 


¥3 

Ri 
Cn Ry 
| CD] [Rs 
Ay iT 


Figure 6.4: Multi-compartment electrical model of a rendered dendrite fiber. 


Our multi-compartment model reveals a 3 branch, 10 node, 27 edge structure to the fiber. Note that 
we have included the Nernst potentials (Definition: "Nernst potentials", p. 12), the nervous impulse as a 
current source, and the additional leftmost edges depicting stimulus current shunted by the cell body. 

We will continue using our previous notation, namely: R; and R,, denoting cell body (Definition: "axial 
resistance", p. 6) and membrane (Definition: "membrane resistance", p. 6) resistances, respectively; x 
representing the vector of potentials 71 ...% 19, and y denoting the vector of currents y,...y27. Using the 
typical value for a cell’s membrane capacitance, 


c= (uF /em?), 
we derive (see variable conventions (Section 2.1.1: Nerve Fibers and the Strang Quartet)): 
Definition 6.5: Capacitance of a Single Compartment 


Cy = 2na-e 


This capacitance is modeled in parallel with the cell’s membrane resistance. Additionally, letting Ac, 
denote the cell body’s surface area, we recall that its capacitance and resistance are 


Definition 6.6: Capacitance of cell body 
Cob = Ache 
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Definition 6.7: Resistance of cell body 
Rep = AcbPm- 


6.6.1.3 Applying the Strang Quartet 
6.6.1.3.1 Step (S1’)—Voltage Drops 


Let’s begin filling out the Strang Quartet. For Step (S1’), we first observe the voltage drops in the figure. 
Since there are a whopping 27 of them, we include only the first six, which are slightly more than we need 


to cover all variations in the set: 
€y =X 


eg = 2%, — En, 


€6 = %2— V3... 


€27 = 110 — En 
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In matrix for, letting b denote the vector of batteries, 


e=b—Az where b= (—Em) 


EF OC OO FF O CO FP OO CO FP OO CO rF OO OF OO Cm OlULmr RLU COUUCUOCUlUDmRRECLUCUCcCOUlUCUOUC Sl 
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and 
-1 0 0 0 0 0 0 0 0 
—1 0 0 0 0 0 0 0 0 
-l1 1 0 0 0 0 0 0 0 0 
0 -l1 0 0 0 0 0 0 0 0 
0 -l 0 0 0 0 0 0 0 0 
0 -l 1 0 0 0 0 0 0 0 
0 0 -1 0 0 0 0 0 0 0 
0 O -1 0 0 0 0 0 0 0 
0 O -1 1 0 0 0 0 0 0 
0 0 O 1 -1 0 0 0 0 0 
0 0 O O -1 0 0 0 0 0 
0 0 O O -1 0 0 0 0 0 
0 0 O O -1 1 0 0 0 0 
A= 0 0 O O 0 -1 0 0 0 0 
0 0 0 O 0 -1 0 0 0 0 
0 0 O O 0 -1 1 0 O 0 
0 0 O 0 0 0 -1 0 0 0 
0 0 0 0 0 0 -1 0 0 0 
0 0 O -1 0 0 0 +1 +0 0 
0 0 0 0 0 0 0 -1 0 0 
0 0 0 0 0 0 0 -1 0 0 
0 0 0 0 0 0 0 -1 1 = 0 
0 0 0 0 0 0 0 0 -1 0 
0 0 0 0 0 0 0 0 -1 0 
0 0 0 0 0 0 0 0 -1 #1 
0 0 0 0 0 0 0 0 0 =! 
0 0 0 0 0 0 0 0 0 =! 


Although our adjacency matrix A is appreciably larger than our previous examples, we have captured the 
same phenomena as before. 


6.6.1.3.2 Applying (S2): Ohm’s Law Augmented with Voltage-Current Law for Capacitors 


Now, recalling Ohm’s Law and remembering that the current through a capacitor varies proportionately 
with the time rate of change of the potential across it, we assemble our vector of currents. As before, we list 
only enough of the 27 currents to fully characterize the set: 


de, 
Oe 
Y1 cb dt 


€2 


~ Reb 


Y2 
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In matrix terms, this compiles to 


where 


69 


Conductance matrix 


Co Oo O&O 
Co Oo O&O 
Co Oo O&O 
Co Oo O&O 
Co Oo Oo 
Co Oo Oo 
Co Oo O&O 
Co Oo Oo 
Co Oo O&O 
Co Oo Oo 
Co Oo O&O 
co Oo O&O 
Co Oo O&O 
SS oo oS 
Co Oo O&O 
Co Oo Oo 
Co Oo Oo 
Co Oo O&O 
Co Oo O&O 
Co Oo Oo 


and 
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Capacitance matrix 


Cag 


oS 


SS 3S Sa SOF OS oS Sa Sa a Ss OS OS: a EOS oO TOS. POS OO, (OS OS Oar OS 
SS Oar Sat ae SQ Se © SO SE MS SS 8S Oar 2S IS FS) 8a wa Er IS SS SS SO: OS SS 
S BSS e'S: 6 eS Se ose) SS SB SS eS. Se) See es a eee 


Se Ce Be | SS Ss SI SS CB BS SS SS  .O VO O  O OS es OS 


ce > a > Te CDS = | ES DS ED ED | SED | a = SED I CE > DO SD > CD OD =) 


Or eS eS Oa OO HO Oa OS SS OS OO OS OD sO Or (Oa OS SS Oa SS Ss OO eS 
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OS OO. OOo oS oS aa OO Oe SOS 
[co a a > a ag SD a A | Sd a SD a 
SS Or 1S Oo "OO. Se 3S 1S OS SS. OS 


(0.19)0 


Ce 1) SET ICS OO SS en Sr Er MS sa FS > KOS 


3 
ane oS ae > io A Ca Sle Rs RR SEP OS RC a ED OR OS Ee OS SP 


So eS 6. So SS. e SS 2S oS So 6:o 6 Sf Se, Ss ae Ss 
3 


Cr Or SS" Oa Oa ah a a Oe A SS a HO OS OS 8S OS SO HO SS FS Oe OS OS a SS 
OS SS 1S ES OS MOS OO RS SS OS OO eS "SS OS OS SS OS OS 1 sO OS OS 
SS SD SS CS Se OS a SS OS SD OS SO Ss SS oO OS SE CS SS Ss - OS 
SSS Ser). SS. se: Se Ss Se SSeS Se: SS SSeS 
SP 8 SS Ss Oe I SS OS SO OS SS OO IO OS Oa OO CS OO 
PG S| Br SOs OE SIG Se, 2 Ss SO SS SIS > a IO SSS OS FS OO OS SS. CO 
S.6°o..o 2S eS Seo eos | 

St GQ SS ar SOs. (Sp OS GS Sa eS ED 1 SS. SGD >! I IO. sD CS eS SS DO St C6 
SS) Sy a SS Ss SP a SS IS SS. SC EPs) SS Sy TO eS SS IE. SS FS BS Ss OS 
SS Sa SS BS SS OE SS RE ar Si Sa Se sO. OS CS SP ES ED CS SS > SS CS 
SS: SS Ca SO Sa Os SS SS HO AC BO SE SS SS SS $B OS, SS NOs SS 
Sy st Sa, SS SS SOS Ss "SE RE SS, SC CS S$ SS, SS, CS OS Oe SS > SO HS IOS OS 


Sr SO Os ISS OS SS OI Se 2S 
(com Sas i > SD ee a DS Dy < 
Se fe ei 0S eee: See eS 


SS. “Or oS Sar (a 


Oo OS 1S Oe Se OS Sa Ca Sa SO Sa 2 OS SS HO Se Oe SES OS 1S OO LO OS. OS 


SS CS OS Ca 8 Ss CS SS St Oe CS Oa Sar 3S OS SS OS OS Ot OS 


SS aS Qe S See Se eS SS a SS Se SSS SS 
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6.6.1.3.3 Step (S3): Applying Kirchoff’s Law 


Our next step is to write out the equations for Kirchoff’s Current Law. We see: 
MW yy S0 
¥3 — ¥4— Y5 — Yo = 0 
Yo — ¥7 — Ys — Yo = 9 

Yo — Yio — Yio = 0 
Yo — Y11 — Yr2 — Yi3s = 0 
Yo — Y11 — Y1r2 — iz = 0 
Ys — Ya — Yi5 — Yio =O 

Y16 — Yi7 — Yis = 0 
Yio — Y20 — Yar — yo2 = 0 
Y22 — Y23 — Yaa — yos = 0 

Y25 — Y26 — Y27 = 0 


Since the B coefficient matrix we’d form here is equal to A’, we can say in matrix terms: 


where the vector f is composed of f; = %9 and fo...27 = 0. 


6.6.1.3.4 Step (S4): Stirring the Ingredients Together 


Step (S4) directs us to assemble our previous toils together into a final equation, which we will then endeavor 
to solve. Using the process (Section 6.1.2.2.4: Step (S4): Assembling) derived in the dynamic Strang module, 
we arrive at the equation 


b 

AToae + ATGAx = ATGD+ f+ aoe, 

which is the general form for RC circuit potential equations. As we have mentioned, this equation presumes 
knowledge of the initial value of each of the potentials, x (0) = X. 


(6.20) 
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Observing our circuit (Figure 6.4: 3-branch Dendrite Model), and letting Ee = Gioo, we calculate the 
13 ) 


necessary quantities to fill out (6.20)’s pieces (for these calculations, see dendrite.m 


Cg 0° 20, Oy 08 se 408. 20>. 30! 0 
0: =e OO, OC Ooh Oe. 40 
O° 80: 26. OS 0? Oi OR OE 
O <<. O05 0 Oo -M- aOs° O> 70 
Pore O° 0) OO O-O% 16 OF Or “O° -O 
GO. “Ge 40... CE: O- OA 
OF <0). “HOS 0: OP sO - Oe 0 
0 0 0 0 0 0 0 Gr 0 0 
0 20: 90-10) Q)- Gy. i > “OR: 0 
Oe. 9 0. - GL SOS GR, 405° ° Shu SO i. 
Gi+Gn -Gi 0 0 0 0 0 0 
—G; 2G; +Gm —G; 0 0 0 0 0 
0 —-G; 2G;+Gn —G; 0 0 0 0 
0 0 —G; 3G; —G; 0 0 —G; 
Aap 0 0 0 —G; 2G;+Gn —-G; 0 0 
0 0 0 —G; 2G;+ Gn —G; 0 
0 0 0 0 —G; G,+Gr, 0 
0 0 0 —G; 0 0 0 2G; + Gn 
0 0 0 0 0 0 —G; 
0 0 0 0 0 0 0 
Geb 
Gm 
Gm 
0 
ATGb = Em Gm 
Gm 
Gm 
Gm 
Gm 
Gin 


'3http://www.ece.rice.edu/~rainking/dendrite/matlab/dendrite.m 
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and an initial (rest) potential of 


xz (0) = Em 


[a a CC oe oe 


6.6.1.4 Applying the Backward-Euler Method 


Since our system is so large, the Backward-Euler method is the best path to a solution. Looking at the 
matrix ATCA, we observe that it is singular and therefore non-invertible. This singularity arises from the 
node connecting the three branches of the fiber and prevents us from using the simple equation x’ = Bx+g, 
we used in earlier Backward-Euler-ings (Section 6.4). However, we will see that a modest generalization to 


our previous form yields (6.21): 
De' =Exr+g (6.21) 


capturing the form of our system and allowing us to solve for x(t). We manipulate (6.21) as follows: 
Da’ =Ex+q 


& (t) — &(t — dt) 


- dt 


= Fz(t)+g 
(D — Edt)  (t) = Dz (t — dt) + gdt 
& (t) = (D — Edt)‘ (D& (t — dt) + gdt), 
where in our case 
D=A'CA, 
E =— (A’GA), and 
g=A™Gb+f. 


This method is implemented in dendrite.m'* with typical cell dimensions and resistivity properties, yielding 
the following graph of potentials. 


'4http://it.is.rice.edu/~rainking /dendrite/matlab/dendrite.m 
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Figure 6.5: 


Graph of Dendrite Potentials 
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(a) Large view of potentials. (b) Zoomed view of potentials showing maxima. 


Chapter 7 


Complex Analysis 1 


7.1 Complex Numbers, Vectors and Matrices’ 


7.1.1 Complex Numbers 


A complex number is simply a pair of real numbers. In order to stress however that the two arithmetics 
differ we separate the two real pieces by the symbol 7. More precisely, each complex number, z, may be 
uniquely expressed by the combination «+ 7zy, where x and y are real and 7 denotes \/—1. We call «x the real 


part and y the imaginary part of z. We now summarize the main rules of complex arithmetic. 
If z] = 2, + 2y, and z2 = x2 + iy2 then 
Definition 7.1: Complex Addition 
a+ 22 = 4, + 22 +i(y1 + ye) 
Definition 7.2: Complex Multiplication 
222 = (%1 + tyr) (@2 + ty2) = 21%2 — yrye + 4 (w1y2 + L2y1) 
Definition 7.3: Complex Conjugation 
A =r —- Wy 
Definition 7.4: Complex Division 


Zz. — 2129 _ £1%2+y1y2+i(r2yi—21y2) 


22 0 Ba Ba 227 +27 


Definition 7.5: Magnitude of a Complex Number 


laa) = Vaz = fai? t+ yn? 


7.1.2 Polar Representation 
In addition to the Cartesian representation z = x + iy one also has the polar form 
z = |z| (cos (@) + isin (6)) 


where @ = arctan (4). 
This form is especially convenient with regards to multiplication. More precisely, 


2122 = |21\|zZ2| (cos (@,) cos (G2) — sin (61) sin (62) + ¢ (cos (61) sin (02) + sin (6;) cos (@2))) 
|z1||zZ2| (cos (0, + 02) + isin (01 + 02)) 


As a result: 
2” = (|z|)” (cos (n) + isin (né)) 


1This content is available online at <http://cnx.org/content /m10504/2.5/>. 


75 


(7.1) 


76 CHAPTER 7. COMPLEX ANALYSIS 1 


7.1.3 Complex Vectors and Matrices 


A complex vector (matrix) is simply a vector (matrix) of complex numbers. Vector and matrix addition 
proceed, as in the real case, from elementwise addition. The dot or inner product of two complex vectors 
requires, however, a little modification. This is evident when we try to use the old notion to define the length 
of a complex vector. To wit, note that. if: 


then 


z?z=(1+i)? + (1-4)? =14+2i-141-21-1=0 


Now length should measure the distance from a point to the origin and should only be zero for the zero 
vector. The fix, as you have probably guessed, is to sum the squares of the magnitudes of the components 
of z. This is accomplished by simply conjugating one of the vectors. Namely, we define the length of a 
complex vector via: 


(eV 2 (7.4) 


In the example above this produces 


Vata? + (La)? = va =2 


As each real number is the conjugate of itself, this new definition subsumes its real counterpart. 
The notion of magnitude also gives us a way to define limits and hence will permit us to introduce 


1 


complex calculus. We say that the sequence of complex numbers, ¢ 2, | n = 2 , converges to the 


complex number zo and write 
Zn — 20 


or 


zy = limit z,, 
n— co 


when, presented with any € > 0 one can produce an integer N for which |Zn — zo] < € when n > N. As an 
example, we note that (2)" — 0. 


7.1.4 Examples 


Example 7.1 
As an example both of a complex matrix and some of the rules of complex arithmetic, let us 
examine the following matrix: 


(7.5) 


zy 
II 
Se Be Se 


Let us attempt to find FF. One option is simply to multiply the two matrices by brute force, 
but this particular matrix has some remarkable qualities that make the job significantly easier. 
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Specifically, we can note that every element not on the diagonal of the resultant matrix is equal to 
0. Furthermore, each element on the diagonal is 4. Hence, we quickly arrive at the matrix 


(7.6) 


| 
oO Oo UF 


0 
0 
4 
0 


nn) 


0 

4 

0 

0 
= Aj 

This final observation, that this matrix multiplied by its transpose yields a constant times the 

identity matrix, is indeed remarkable. This particular matrix is an example of a Fourier matrix, 


and enjoys a number of interesting properties. The property outlined above can be generalized for 
any F,,, where F refers to a Fourier matrix with n rows and columns: 


F, FPF, = nl (7.7) 


7.2 Complex Functions’ 


7.2.1 Complex Functions 


A complex function is merely a rule for assigning certain complex numbers to other complex numbers. The 
simplest (nonconstant) assignment is the identity function f(z) = z. Perhaps the next simplest function 
assigns to each number its square, i.e., f(z) = z?. As we decomposed the argument of f, namely z, into 
its real and imaginary parts, we shall also find it convenient to partition the value of f, z? in this case, into 
its real and imaginary parts. In general, we write 


fw +ty) = u(a,y) + iv (zy) (7.8) 
where u and v are both real-valued functions of two real variables. In the case that f (z) = 2? we find 
u(e,y)=a°—y? 


and 
v (x,y) = 2ay 


With the tools of complex numbers (Section 7.1), we may produce complex polynomials 
f(g=2™ + ene te Faz teg (7.9) 


We say that such an f is order m. We shall often find it convenient to represent polynomials as the product 


of their factors, namely 
f@=6-—AD™@— 2)? eA)” (7.10) 


Each A; is a root of degree d;. Here h is the number of distinct roots of f. We call \; a simple root when 
d; = 1. We can observe the appearance of ratios of polynomials or so called rational functions. Suppose 


?This content is available online at <http://cnx.org/content/m10505/2.7/>. 
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in rational, that f is of order at most m— 1 while g is of order m with the simple roots {A1,...,Am}. It 
should come as no surprise that such gq should admit a Partial Fraction Expansion 


a(z)=>> Ear (7.11) 


One uncovers the q; by first multiplying each side by z— A; and then letting z tend to \,;. For example, if 


1 71 q2 
= 7.12 
z+1  z41 a z—-1 ( ) 
then multiplying each side by z + 7 produces 
l : 
see lea (7.13) 
z—4 z—1 
Now, in order to isolate q, it is clear that we should set z = —7. So doing we find that q, = 4 In order to 
find q2 we multiply (7.12) by z —7. and then set z = 7. So doing we find qz = 3, and so 
13 z 
—_ seh : (7.14) 
a a z+ zZ—4 
. Returning to the general case, we encode the above in the simple formula 
qj = limit (z—A;) a (z) (7.15) 
You should be able to use this to confirm that 
1/2 1/2 
ete cae (7.16) 


z+1 oz+ti 2-1 


Recall that the transfer function we met in The Laplace Transform (p. 58) module was in fact a matrix 
of rational functions. Now, the partial fraction expansion of a matrix of rational functions is simply the 
matrix of partial fraction expansions of each of its elements. This is easier done than said. For example, the 
transfer function of 


0 1 
B= 
—-l1 0 
is 
-1 z 1 
(zI — B) = i 
—l z 
; : (7.17) 
if i) i (ie F 
PEON ghee. HD ne as 


The first line comes form either Gauss-Jordan by hand or via the symbolic toolbox in Matlab. More 
importantly, the second line is simply an amalgamation of (7.12) and (7.14). Complex matrices have finally 
entered the picture. We shall devote all of Chapter 10 to uncovering the remarkable properties enjoyed 
by the matrices that appear in the partial fraction expansion of (zJ — B)" Have you noticed that, in our 
example, the two matrices are each projections, and they sum to J, and that their product is 0? Could this 
be an accident? 

In The Laplace Transform (Section 6.2) module we were confronted with the complex exponential. By 
analogy to the real exponential we define 
n 


on = z 
v=)> at (7.18) 
n=0 
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and find that 


ee pag Oe GOs ek Oe a. 
2 4 ; 3 5 
= 1-$4+% +i(8 e+ ey, (7.19) 


= cos(@) + isin (8) 


With this observation, the polar form (Section 7.1.2: Polar Representation) is now simply z = |z\e’?. 
One may just as easily verify that 
ci 4 (a0 


cos (0) = 5 


and ; 
ef = e(- 8 


sin (0) = ry 


These suggest the definitions, for complex z, of 


iz (i-—z) 
a ei (7.20) 


cos (z) 
and ; 
ei? — e(-H)z 


i (7.21) 


As in the real case the exponential enjoys the property that 
e71+22 — e%1¢%2 


and in particular 
erry = ere 


(7.22) 


I 


e*cos (y) + ie*sin (y) 


Finally, the inverse of the complex exponential is the complex logarithm, 
In (z) = In({z|) +76 (7.23) 


for z = |z|e*®. One finds that In (—1 + i) = In (V2) + i. 


7.3 Complex Differentiation® 


7.3.1 Complex Differentiation 


The complex f is said to be differentiable at zo if 


rt £—F ) 


Z—Zo ZS ZO 


exists, by which we mean that, 
Ff (én) ~ f (20) 
en — 20 
converges to the same value for every sequence {z,,} that converges to zo. In this case we naturally call 
the limit 4 f (zo) 


’This content is available online at <http://cnx.org/content /m10276/2.7/>. 
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To illustrate the concept of ’for every’ mentioned above, we utilize the following picture. We assume 
the point zo is differentiable, which means that any conceivable sequence is going to converge to zo. We 
outline three sequences in the picture: real numbers, imaginary numbers, and a spiral pattern of both. 


Sequences Approaching A Point In The Complex Plane 


Figure 7.1: The green is real, the blue is imaginary, and the red is the spiral. 
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7.3.2 Examples 


Example 7.2 
The derivative of z? is 2z. 


27-297 


limit = limit @-#0)let20) 
Zz «70 ZZ as #0. (7.24) 
= 220 
Example 7.3 
The exponential is its own derivative. 
limit <=" = limit @—°=1 
zg = 0 z—zq 0 70 
= sie co §6(z-20)" 
= e*limit bie (ati) (7.25) 
= e70 


Example 7.4 
The real part of z is not a differentiable function of z. 

We show that the limit depends on the angle of approach. First, when z, — zo on a line parallel 
to the real axis, e.g., 2, = ®o + + + iyo, we find 


tices 
Os ae) =i (7.26) 


limit a - 
n-co Lg + 7 + tY0 — Lo + 10 


while if z, — zo in the imaginary direction, e.g., 2, = xo +12 (Yo + +), then 


to — Xo 


limit - —=0 7.27 
n=co zy +i (yo +) — ao + iyo ee 


7.3.3 Conclusion 
This last example suggests that when f is differentiable a simple relationship must bind its partial derivatives 
in x and y. 
Proposition 7.1: Partial Derivative Relationship 
If f is differentiable at zo then € f (zo) = 242) = (i242) 
Proof: 
With z = x + iyo, 


ff (0) = limiy Aap 
= limit fei) Tea it) a5 

L— 2X0 

— Of) 

= Ox 


Alternatively, when z = 2 + iy then 
f(z)—f (0) 


Z—ZOQ 


I 


tf (zo) limit 
Z— Zo 

— fim Lleottv)—f(totiyo) 
Toni c= (7.29) 
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7.3.4 Cauchy-Reimann Equations 


In terms of the real and imaginary parts of f this result brings the Cauchy-Riemann equations. 


Ou Ov 
and 5 3 
v U 
pes nal eS 31 
Ox Oy eh) 


Regarding the converse proposition we note that when f has continuous partial derivatives in region obeying 
the Cauchy-Reimann equations then f is in fact differentiable in the region. 

We remark that with no more energy than that expended on their real cousins one may uncover the rules 
for differentiating complex sums, products, quotients, and compositions. 

As one important application of the derivative let us attempt to expand in partial fractions a rational 
function whose denominator has a root with degree larger than one. As a warm-up let us try to find q11 
and qi,2 in the expression 

Z+2 My 1,2 
Cal? erie 


Arguing as above, it seems wise to multiply through by (z + 1)? and so arrive at 
z2+2=q1(2+1)4+Q,2 (7.32) 


On setting z = —1 this gives q).2 = 1. With q:,2 computed, (7.32) takes the simple form z+1=q1 (z+ 1) 


and so qi,2 = 1 as well. Hence, 
z+2 1 1 


(+1)? 2+1 (241)? 


This latter step grows more cumbersome for roots of higher degrees. Let us consider 


(+2)? _ Mi, 41,2 4 11,3 


(eh SBt eae? ea )® 


The first step is still correct: multiply through by the factor at its highest degree, here 3. This leaves us 
with 
(+2) =qa(z+l1)P+a2(zt+)tags (7.33) 


Setting z = —1 again produces the last coefficient, here q;,3 = 1. We are left however with one equation in 
two unknowns. Well, not really one equation, for (7.33) is to hold for all z. We exploit this by taking two 
derivatives, with respect to z, of (7.33). This produces 


2(2+2) = 2a 1 (z+1) +412 


and 2 = qi The latter of course needs no comment. We derive qi,2 from the former by setting z = —1. This 
example will permit us to derive a simple expression for the partial fraction expansion of the general 
proper rational function q = f where g has h distinct roots {A,,...,An,} of respective degrees {d,,...,dy}. 
We write 


h 
W@=> > ea (7.34) 


and note, as above, that q;,, is the coefficient of (z — d;)a—* in the ration 


ry (2) =a (2) (z—Ay)” 
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al function 


Hence, q;,, may be computed by setting z = A, in the ratio of the d; — kth derivative of r; to (d; — k)!. 


That is, 
a eee ae {(z yaa} (7.35) 
zt zr; (d; = k)! dzti-k 4 
As a second example, let us take 
1 0 0 
B=] 1 0 (7.36) 
01 1 
and compute the ®,, matrices in the expansion 
ay 0 0 
-1 1 1 1 Ai 
(zI-B) = (x1)(=3) 3 0 = yap + Ge yPi2+ ya 
1 
(2-1)?(z-3) (2-1)(2-3) sz 
The only challenging term is the (3,1) element. We write 
1 on a1 | 1,2 _ 92,1 
G17 G3) 2k @aiy 2-3 
It follows from (7.35) that 
as fa 
a1 = a (251) (7.37) 
= -1/4 
and 
a2 = x3! (7.38) 
= -1/4 
and 
1 \2 
dake (>) 2 (7.39) 
= 1/4 
It now follows that 
; 1 0 0 0 0 0 ‘ 0 0 0 
(zl — B)" =— ]} -1/2 0 0 | +——>3 0 0 f+ 1/2 1 O (7.40) 
z—1 (z-1) z—3 
—1/4 -1/2 1 -1/2 0 0 1/4 1/2 0 


In closing, let us remark that the method of partial fraction expansions has been implemented in Matlab. 


In fact, (7.37), (7.38), and (7.39) all follow from the single command: [r,p, 
7 -3]). The first input argument is Matlab-speak for the polynomial f (z) 
corresponds to the denominator 


g(z) = (2-1)? (2-3) = A -52? 472-3 


k]=residue([0O 0 0 1],[1 -5 
= 1 while the second argument 
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7.4 Exercises: Complex Numbers, Vectors, and Functions’ 


7.4.1 Exercises 


Exercise 7.1 (Solution on p. 85.) 
Express |e*| in terms of x and/or y. 

Exercise 7.2 (Solution on p. 85.) 
Confirm that e'™) = z and In (e”) = z. 

Exercise 7.3 (Solution on p. 85.) 


Find the real and imaginary parts of cos(z) and sin (z). Express your answers in terms of regular 
and hyperbolic trigonometric functions. 


Exercise 7.4 (Solution on p. 85.) 
Show that cos? (z) + sin? (z) = 1 
Exercise 7.5 (Solution on p. 85.) 


With z” = e”™©) for complex z and w compute Vi 

Exercise 7.6 (Solution on p. 85.) 
Verify that cos(z) and sin(z) satisfy the Cauchy-Riemann equations and use the proposition to 
evaluate their derivatives. 

Exercise 7.7 (Solution on p. 85.) 


Submit a Matlab diary documenting your use of residue in the partial fraction expansion of the 
transfer function of 


2 0 O 
B=|-1 4 0O 
O -1 2 


4This content is available online at <http://cnx.org/content /m10506/2.5/>. 


Solutions to Exercises in Chapter 7 


Solution to Exercise 7.1 (p. 84) 
Pending completion of assignment. 
Solution to Exercise 7.2 (p. 84) 
Pending completion of assignment. 
Solution to Exercise 7.3 (p. 84) 
Pending completion of assignment. 
Solution to Exercise 7.4 (p. 84) 
Pending completion of assignment. 
Solution to Exercise 7.5 (p. 84) 
Pending completion of assignment. 
Solution to Exercise 7.6 (p. 84) 
Pending completion of assignment. 
Solution to Exercise 7.7 (p. 84) 
Pending completion of assignment. 
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Chapter 8 


Complex Analysis 2 


8.1 Cauchy’s Theorem’ 


8.1.1 Introduction 


Our main goal is a better understanding of the partial fraction expansion of a given transfer function. With 
respect to the example that closed the discussion of complex differentiation, see this equation (7.36) - In this 


equation (7.40), we found 


a 1 1 
I-B)'= P, D P. 
CS ee ig seye ot za 
where the P; and D; enjoy the amazing properties 
1. 
BP, = PRB 
: ; (8.1) 
= AP +D, 
and 
BP, = P2B = d2P» 
2. 
P+P=I (8.2) 
P? =P, 
Pp? = Py 
and 
D,? =0 
3: 
PD, = DP. 
1Yi if (8.3) 
= Di, 
and 


P2,D, = Di P, = 0 


In order to show that this always happens, i.e., that it is not a quirk produced by the particular B in this 
equation (7.36), we require a few additional tools from the theory of complex variables. In particular, we 
need the fact that partial fraction expansions may be carried out through complex integration. 


!This content is available online at <http://cnx.org/content /m10264/2.8/>. 
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8.1.2 Integration of Complex Functions Over Complex Curves 


We shall be integrating complex functions over complex curves. Such a curve is parameterized by one 
complex valued or, equivalently, two real valued, function(s) of a real parameter (typically denoted by ¢). 
More precisely, 

C={2@) =e) +y@)|)a<t< 6} 


For example, if x(t) = y(t) =t while a = 0 and b = 1, then C is the line segment joining 0+ 70 to 1 +7. 


We now define 
b 
fi@usf f(2(t)) 2! (t) at 


For example, if C = {t+ it |0<t< 1} as above and f(z) =z then 


foee fi (+ it) (1 +1) dt oP t—t+i2tdt =% 


while if C is the unit circle {e | 0 <t < 27} then 


20 20 20 
: zdz = ‘ ete dt = if Bete = if cos (2t) + isin (2t) dt = 0 
0 0 0 


Remaining with the unit circle but now integrating f (z) = + we find 


20 
peta = e ietdt = Qri 
0 


We generalize this calculation to arbitrary (integer) powers over arbitrary circles. More precisely, for 
integer m and fixed complex a we integrate (z — a)" over 


C(a,r) = {at+re*|0<t<2nr} 
the circle of radius r centered at a. We find 


f (z-—a)"dz = fig (a + ret — a)" rie''dt 


oe (8.4) 
= gpmtl ie et(m+1)t qe 


2ri if m=—-l 


, (z—a)"dz = ir i cos ((m + 1) t) + isin ((m + 1) t) dt = 
0 


0 otherwise 


When integrating more general functions it is often convenient to express the integral in terms of its real 
and imaginary parts. More precisely 


[f@ae= ful) +ivende +i f u(ey) + iv(ey dy 


[i@a=fucuac- f o(eydy+i f v(eyacti f uleydy 


J f@dz = fu@@,y@)a() - v(e(t),y@)y Wat + 
ifpu(a(t),y@)y (Wt) +(e (t),y@)) 2" (© dt 


The second line should invoke memories of: 
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Theorem 8.1: Green’s Theorem 
If C is a closed curve and M and N are continuously differentiable real-valued functions on Cin, 
the region enclosed by C, then 


[Mace fxay= [ eG 


Applying this to the situation above, we find, so long as C is closed, that 


[i@a=- [= Seedy +i ff 45 5y ty 


At first glance it appears that Green’s Theorem only serves to muddy the waters. Recalling the Cauchy- 
Riemann equations (Section 7.3.4: Cauchy-Reimann Equations) however we find that each of these double 
integrals is in fact identically zero! In brief, we have proven: 


Theorem 8.2: Cauchy’s Theorem 
If f is differentiable on and in the closed curve C then f{ f (z) dz =0. 


Strictly speaking, in order to invoke Green’s Theorem we require not only that f be differentiable but 
that its derivative in fact be continuous. This however is simply a limitation of our simple mode of proof; 
Cauchy’s Theorem is true as stated. 

This theorem, together with (8.4), permits us to integrate every proper rational function. More precisely, 
ifq= Pa where f is a polynomial of degree at most m—1 and g is an mth degree polynomial with h distinct 
zeros at {A; | j = {1,...,h}} with respective multiplicities of {m, | 7 ={1,...,h}} we found that 


> re eee luaees (8.5) 
j=l k= i ( — Aj) 


Observe now that if we choose r; so small that ,; is the only zero of g encircled by C; = C (Aj, 1r;) then by 


Cauchy’s Theorem 
mi 1 
q(z)dz= y wn | dz 
i k=1 (z = rF 


In (8.4) we found that each, save the first, of the integrals under the sum is in fact zero. Hence, 


fa (z) dz = 2nig;1 (8.6) 


With q;,1 in hand, say from this equation (7.35) or residue, one may view (8.6) as a means for computing 
the indicated integral. The opposite reading, i.e., that the integral is a convenient means of expressing qj,1, 
will prove just as useful. With that in mind, we note that the remaining residues may be computed as 
integrals of the product of g and the appropriate factor. More precisely, 


Ja (z) (z- dj)* "dz = 2nig;,k (8.7) 
One may be led to believe that the precision of this result is due to the very special choice of curve and 
function. We shall see ... 


8.2 Cauchy’s Integral Formula’ 


8.2.1 The Residue Theorem 


After Cauchy’s Theorem eqn6 (8.6) and Cauchy’s Theorem eqn7 (8.7) perhaps the most useful consequence 
of Cauchy’s Theorem (Theorem 8.2, Cauchy’s Theorem, p. 89) is the 


?This content is available online at <http://cnx.org/content /m10246/2.8/>. 
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Lemma 8.1: The Curve Replacement Lemma 
Suppose that C2 is a closed curve that lies inside the region encircled by the closed curve C\. If f 
is differentiable in the annular region outside C2 and inside C, then 


[i@a=f rea 
Proof: 


With reference to Figure 8.1 (Curve Replacement Figure) we introduce two vertical segments and 

define the closed curves C3 = abcda (where the be arc is clockwise and the da arc is counter- 
clockwise) and C4 = adcba (where the ad arc is counter-clockwise and the cb arc is clockwise). By 
merely following the arrows we learn that 


[i@a=f reat froat [toa 


As Cauchy’s Theorem (Theorem 8.2, Cauchy’s Theorem, p. 89) implies that the integrals over C3 
and C4 each vanish, we have our result. 


Curve Replacement Figure 


Figure 8.1: The Curve Replacement Lemma 


91 


This Lemma says that in order to integrate a function it suffices to integrate it over regions where it is 
singular, i.e. nondifferentiable. 
Let us apply this reasoning to the integral 


| ee” 


where C encircles both A, and Az as depicted in Figure 8.2. We find that 


bore Seen ee 


Developing the integrand in partial fractions we find 


/ z F) = Al i, 1 d ; AQ - 1 a i QTiAY 
(z — Ax) (2 — Az) at aa Fes Oey z—X»o am ae 


/ z d Qrir2 
ia 
(z — A1) (z — A2) Ag — Ay 


Putting things back together we find 


Similarly, 


z : r A 
if = eam = ani (5 Ab + a) 


= 2ni 


(8.8) 


Image not finished 


Figure 8.2: Concentrating on the poles. 


We may view (8.8) as a special instance of integrating a rational function around a curve that encircles 
all of the zeros of its denominator. In particular, recalling that cauchy’s theorem eqn5 (8.5) and Cauchy’s 
theorem eqn6 (8.6), we find 


h mi jk 
Ja(zjdz = yo a Goi dz 


(8.9) 
= Qt ae qj,1 


To take a slightly more complicated example let us integrate i) over some closed curve C’ inside of 


which f is differentiable and a resides. Our Curve Replacement Lemma now permits us to claim that 


za 


It appears that one can go no further without specifying f. The alert reader however recognizes that in the 
integral over C'(a,r) is independent of r and so proceeds to let r — 0, in which case z — a and f (z) — f (a). 
Computing the integral of along the way we are led to the hope that 


zZ—-a 


LA) 5, = f (a) 2m 


za 
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In support of this conclusion we note that 
z a 1 = 
Lien f HO HO Liha yo f our { aL 


Now the first term is f (a) 277 regardless of r while, as r — 0, the integrand of the second term approaches 
if (a) and the region of integration approaches the point a. Regarding this second term, as the integrand 
remains bounded as the perimeter of C'(a,r) approaches zero the value of the integral must itself be zero. 
This result if typically known as 


Rule 8.1: Cauchy’s Integral Formula 
If f is differentiable on and in the closed curve C' then 


r(a)= 5; f Pa: 


277 zZ-a 


(8.10) 


for each a lying inside C. 


The consequences of such a formula run far and deep. We shall delve into only one or two. First, we note 
that, as a does not lie on C, the right hand side is a perfectly smooth function of a. Hence, differentiating 


each side, we find 
d 1 f(z) 
dat (2) 2% / (z—a)? os can 


for each a lying inside C. Applying this reasoning n times we arrive at a formula for the n-th derivative of 
f at a, 


d” n! z 
dart) = 9 | ae ae 


Td 


for each a lying inside C. The upshot is that once f is shown to be differentiable it must in fact be infinitely 
differentiable. As a simple extension let us consider 


1 ff 
Qni i GIs we 


where f is still assumed differentiable on and in C' and that C encircles both A; and Az. By the curve 
replacement lemma this integral is the sum 


1 f (2) 1 (z) 
mi iy (z— 1) (z- me” * Oni ; (z—A)(z- me 
f(z) 
z—X2 
1 | f (2) jpee Tt OW) 
Qmi J (z— 1) (z— rz)” (Ai — Az)? 


where \; now lies in only C;. As is well behaved in C, we may use (8.10) to conclude that 


ie) is well behaved in Cy we may use (8.11) to conclude that 


Similarly, as ee 
1 f (z) 2 d ( f (a) ) 
ant / (z p=. A1) (z — tae ~ da a- AL here 


These calculations can be read as a concrete instance of 
Theorem 8.3: The Residue Theorem 
If g is a polynomial with roots {A; | j = {1,...,h}} of degree {d; | 7 ={1,...,h}} and Cis a 
closed curve encircling each of the A; and f is differentiable on and in C' then 


h 


INA) Tt re ; 
[ige= Sores (Ay) 


j=1 
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where 


— d; "(z 
1 d4#~*(z— dj) ie) 


tes i) oe d2eG—t 


is called the residue of f at ;. 


One of the most important instances of this theorem is the formula for the Inverse Laplace Transform. 


8.3 The Inverse Laplace Transform: Complex Integration*® 


8.3.1 The Inverse Laplace Transform 


If q is a rational function with poles {A, | 7 = {1,...,4}} then the inverse Laplace transform of q is 


~ Qi 


c@i== / q(2) edz (8.13) 


where C is a curve that encloses each of the poles of g. As a result 


£>* (q) (t) = Sores (Ay) (8.14) 


j=l 


Let us put this lovely formula to the test. We take our examples from discussion of the Laplace Transform 
(Section 6.2) and the inverse Laplace Transform (Section 6.3). Let us first compute the inverse Laplace 


Transform of 
gs 
no “sate 


According to (8.14) it is simply the residue of q(z) e** at z = —1, ie., 


res (—1) = limit —te! 


Zl vA 


This closes the circle on the example begun in the discussion of the Laplace Transform (Section 6.2) and 
continued in exercise one for chapter 6. 
For our next example we recall 


0.19 (s? + 1.5 + 0.27) 


L£ (a1 (s)) = 
Le) (s + 1/6)* (s3 + 1.65582 + 0.49788 + 0.0039) 


from the Inverse Laplace Transform (Section 6.3). Using numde, sym2poly and residue, see fib4.m for 
details, returns 


0.0029 
262.8394 
—474.1929 
r= —1.0857 
—9.0930 
—0.3326 
211.3507 
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and 
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—1.3565 
—0.2885 
—0.1667 
Pi =| —0.1667 
—0.1667 
—0.1667 
—0.0100 


You will be asked in the exercises to show that this indeed jibes with the 


a} (t) = 211.35em — (0.0554 + 4.54644? + 1.085¢+474.19)es + 


e400 


oF (262.842cosh () wanaesoeae ()) 


achieved in the Laplace Transform (Section 6.2) via ilaplace. 


8.4 


1. 


Exercises: Complex Integration’ 


Let us confirm the representation of this Cauchy’s Theorem equation (8.7) in the matrix case. More 
precisely, if ® (z) = (zI — B)~' is the transfer function associated with B then this Cauchy’s Theorem 
equation (8.7) states that 


where 


®), = — | ®(z)(z—d;)* "dz (8.15) 


Compute the ®;, per (8.15) for the B in this equation (7.36) from the discussion of Complex Dif- 
ferentiation. Confirm that they agree with those appearing in this equation (7.40) from the Complex 
Differentiation discussion. 


. Use this inverse Laplace Transform equation (8.14) to compute the inverse Laplace transform of 


1 
s?+2s+2° 


. Use the result of the previous exercise to solve, via the Laplace transform, the differential equation 


d 


ae (x) (t) +. 2(t)=e7'sin(t), «(0) =0 


Hint: Take the Laplace transform of each side. 


. Explain how one gets from r; and p; to @; (t). 
. Compute, as in fib4.m, the residues of L (x2 (s)) and L (a3 (s)) and confirm that they give rise to the 


x2 (t) and x3 (t) you derived in the discussion of Chapter 1 (Section 2.1). 


4This content is available online at <http://cnx.org/content /m10524/2.4/>. 


Chapter 9 


The Eigenvalue Problem 


9.1 Introduction’ 


9.1.1 Introduction 


Harking back to our previous discussion of The Laplace Transform (Section 6.2) we labeled the complex 
number A an eigenvalue of B if AJ — B was not invertible. In order to find such X one has only to find 
those s for which (sf — B)~' is not defined. To take a concrete example we note that if 


1 1 0 
B=|0 1 0 (9.1) 
0 0 2 
then 
: (s —1)(s— 2) s—2 0 
(sf — B)-* = ——,——_ 0 (s—1)(s—2) 0 (9.2) 
(s — 1)" (s— 2) . aig 


and so A; = 1 and 2 = 2 are the two eigenvalues of B. Now, to say that A,;/ — B is not invertible is to 
say that its columns are linearly dependent, or, equivalently, that the null space -V (A;J — B) contains more 
than just the zero vector. We call YW (A;I — B) the jth eigenspace and call each of its nonzero members a 
jth eigenvector. The dimension of VY (A,;I — B) is referred to as the geometric multiplicity of \;. With 
respect to B above, we compute -V (AJ — B) by solving (J — B) x = 0, ie., 


—-l1 0 Vy 0 
0 0 O LQ =] 0 
0 0 1 x3 0 


Clearly 
a6 
W(t ~B)={e( 1 0 0) cer} 


Arguing along the same lines we also find 


HV (nal ~ B)={ e( 0 0 1) leer} 


!This content is available online at <http://cnx.org/content /m10405/2.4/>. 
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That B is 3x3 but possesses only 2 linearly eigenvectors leads us to speak of B as defective. The cause of 
its defect is most likely the fact that A, is a double pole of (sI — B)** In order to flesh out that remark and 
uncover the missing eigenvector we must take a much closer look at the transfer function 


R(s) =(sI — B)' 


In the mathematical literature this quantity is typically referred to as the Resolvent of B. 


9.2 The Resolvent’ 


9.2.1 The Transfer Function 


One means by which to come to grips with R(s) is to treat it as the matrix analog of the scalar function 


1 
s—b 


(9.3) 


This function is a scaled version of the even simpler function ae This latter function satisfies the identity 


(just multiply across by 1 — z to check it) 


1 2 n—1 2" 
— =l4+z4+27°+..4+2 + 
l-<z l-<z 


(9.4) 


for each positive integer n. Furthermore, if |z| < 1 then z” — 0 as n — oo and so (9.4) becomes, in the 
limit, 


1 <a 
ree oe 
n=0 
the familiar geometric series. Returning to (9.3) we write 
1 4 1 b acre aml 
=§ = +5H.4 
eo Da ae 8 s” s™s—b 


and hence, so long as |s| > |b| we find, 


1 ae aa 
ee, 


This same line of reasoning may be applied in the matrix case. That is, 


(t- By? =o4(1 ee = Load BES —~(sh 8") (9.5) 


and hence, so long as |s| >|| B || where || B || is the magnitude of the largest element of B, we find 


(sf —B)7'=s57! % (2) (9.6) 


Although (9.6) is indeed a formula for the transfer function you may, regarding computation, not find it 
any more attractive than the Gauss-Jordan method. We view (9.6) however as an analytical rather than 
computational tool. More precisely, it facilitates the computation of integrals of R(s). For example, if C, is 
the circle of radius p centered at the origin and p >|| B || then 


Jo, (sf —B)7'ds = Beery oy Jo, g's 


(9.7) 
Q2rit 


I 
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This result is essential to our study of the eigenvalue problem. As are the two resolvent identities. Regarding 


the first we deduce from the simple observation 


(sol — B)~* — (s,I — B)"' = (821 — B)* (s,|I — B— 821 + B) (s:I —- B)™* 


that 
R (s2) — R(s1) = (81 — $2) R(s2) R(s1) 


The second identity is simply a rewriting of 
(sI — B) (sf — B)~' =(sI-— B)"'(sI- B) =I 


namely, 
BR(s) = R(s)B 
= sR(s)- 


9.3 The Partial Fraction Expansion of the Resolvent* 


9.3.1 Partial Fraction Expansion of the Transfer Function 


(9.8) 


(9.9) 


The Gauss-Jordan method informs us that R will be a matrix of rational functions with a common denom- 
inator. In keeping with the notation of the previous chapters, we assume the denominator to have the h 


distinct roots, {A; | j = {1,...,h}} with associated multiplicities {m, | 7 = {1,...,A}}. 
Now, assembling the partial fraction expansions of each element of R we arrive at 


h M5 


So) rsa 


jal kal ( 


where, recalling the equation from Cauchy’s Theorem (8.7), the matrix R,;, equals the following: 


— 1 
Ryp = 3; | RO) (z—X j dz 


Example 9.1: Concrete Example 


(9.10) 


(9.11) 


As we look at this example with respect to the eigenvalue problem eqn1 (9.1) and eqn2 (9.2), we 


1 0 0 0 
0 and = =R21 = 0 0 0 
0 0 0 1 


Co Oo 2 


One notes immediately that these matrices enjoy some amazing properties. For example 


Ry? =Ria, Roa? =Rei, RiaRe1=0, and Ri? =0 


Below we will now show that this is no accident. As a consequence of (9.11) and the first resolvent 


identity, we shall find that these results are true in general. 


Proposition 9.1: 
Ra? = R;,1 as seen above in (9.12). 
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Proof: 
Recall that the C; appearing in (9.11) is any circle about \; that neither touches nor encircles any 
other root. Suppose that C; and C;' are two such circles and C;’ encloses C’;. Now 


1 1 


Ry R (z) dz R (z) dz 


Ry? = ae | ROG | Rewyaw 
Rj? = ae | [ROR wave: 
natn gig f J BBs 


Ra =, (fr f Rcawae— f rw) f a piledw) 


1 
Ry? = af BO dz= Ry 


and so 


We used the first resolvent identity, This Transfer Function eqn (9.8), in moving from the second 
to the third line. In moving from the fourth to the fifth we used only 


1 
7 ——dw = 2ni (9.13) 
w-zZ 


1 
i} dz =0 
w—zZz 


The latter integrates to zero because C; does not encircle w. 
From the definition of orthogonal projections (Definition: "orthogonal projection", p. 51), which 
states that matrices that equal their squares are projections, we adopt the abbreviation 


P; — Rya 


With respect to the product PjP,, for j # k, the calculation runs along the same lines. The 
difference comes in (9.13) where, as C; lies completely outside of Cy, both integrals are zero. 
Hence, 


Proposition 9.2: 
If 7 Ak then Pj P, = 0. 
Along the same lines we define 
D; = Rj2 
and prove 


Proposition 9.3: 

If1<k<m,—1then Dj* = Rypyi. Dj =0. 
Proof: 

For k and | greater than or equal to one, 


ese ere. z) (z—A;)*dz w) (w —d;)'dw 
Rares = oo / R(z)(2-Aj)Rd i R(w) (w — rj) 
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Sears z w) (z — A;)*(w — A;)'dwdz 
RireRyn = ooo ff RO) ROW) dy)*(w — j)!dwa 


Pe eee ae OO eee es 
Ry ri Aj i41 = rae ce (z—A;)"( Aj) dwd 


Ryka = aap f 2 @-r)* f WA tude a pRw w— ay f CoM) joa 
Ry epi Rj = = i R(z)(z— Aj)" dz = Ry egins 
because ; 
/ (WON) ay = 2ni(z — dj)! (9.14) 
and 


k 
[a0 
W— Zz 


With k =1 = 1 we have shown Rj2” = Rj, i.e., Dj” = Rj,3. Similarly, with k = 1 and | = 2 we 
find Rj 2Rj,3 = Rj,4, 1e., D;* = Rj,4. Continuing in this fashion we find Rj ,Rjn41 = Ry x42 = J, 
or Dr = Ry x+2. Finally, at k =m, — 1 this becomes 


be 
Dy? = Rimyt1 = 


! [Re (z — Aj)" dz =0 


Qri 


by Cauchy’s Theorem. 
With this we now have the sought after expansion 


mj-1 


h 
1 
R(z)= P;4 D;* (9.15) 


along with the verification of a number of the properties laid out in Complex Integration eqns 1-3 
(8.1). 


9.4 The Spectral Representation’ 


With just a little bit more work we shall arrive at a similar expansion for B itself. We begin by applying the 
second resolvent identity (9.9) to P;. More precisely, we note that the second resolvent identity (9.9) implies 
that 
BP; = P;B 
(9.16) 
= fo, 2R(2)—Idz 


pp | zR(z) dz 
Cj 


PB= f Ree-Ade+r [Rae 
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Summing this over j we find 
h h h 
BYP) => (dP) + 0D; (9.17) 


We can go one step further, namely the evaluation of the first sum. This stems from the eqn in the discussion 
of the transfer function (9.7) where we integrated R(s) over a circle C, where p >|| B ||. The connection to 
the P; is made by the residue theorem. More precisely, 


Comparing this to the eqn (9.7) from the discussion of the transfer function we find 


h 
bye aes (9.18) 
j=l 


and so (9.17) takes the form 
h h 
B=) °\jP; +5 D; (9.19) 


It is this formula that we refer to as the Spectral Representation of B. To the numerous connections 
between the P; and D; we wish to add one more. We first write (9.16) as 


(BAG) P= D; 
and then raise each side to the m; power. As ie = P; and Dy = 0 we find 
(B— A, 1)" P; =0 (9.20) 


For this reason we call the range of P; the jth generalized eigenspace, call each of its nonzero members 
a jth generalized eigenvector and refer to the dimension of R(P;) as the algebraic multiplicity of 
Aj. With regard to the first example (p. 95) from the discussion of the eigenvalue problem, we note that 
although it has only two linearly independent eigenvectors the span of the associated generalized eigenspaces 
indeed fills out R°. One may view this as a consequence of P, + Py = I, or, perhaps more concretely, as 


T 
appending the generalized first eigenvector ( 01 0 ) to the original two eigenvectors ( 10 0 ) 


and ( 0 0 1 ) . In still other words, the algebraic multiplicities sum to the ambient dimension (here 3), 


while the sum of geometric multiplicities falls short (here 2). 


9.5 The Eigenvalue Problem: Examples’ 


We take a look back at our previous examples in light of the results of two previous sections The Spectral 
Representation (Section 9.4) and The Partial Fraction Expansion of the Transfer Function (Section 9.3). 
With respect to the rotation matrix 


O 1 
-1 0 


B= 
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we recall, see Cauchy’s Theorem eqn6 (8.6), that 


1 S 1 
R(s)= 4 
Cage oe a aoe ee 
1 1/72 = 1 1/2 2 
R(s)= : 2 | +— . z (9.21) 
s—t 4 1/2 sti\ = 1/2 
1 1 
R(s)= Pi+ P. 
(s) Pe Os ee Pa 
and so 
Lo = 1/22 
B=\P,+)2Po =i / 2 | 4 / a. 
i 1/2 3 1/2 


From m, = mg = 1 it follows that R (P,) and R (P2) are actual (as opposed to generalized) eigenspaces. 
These column spaces are easily determined. In particular, R (P,) is the span of 


while R (P2) is the span of 


eg = : 
-i 
To recapitulate, from partial fraction expansion one can read off the projections from which one can read off 
the eigenvectors. The reverse direction, producing projections from eigenvectors, is equally worthwhile. We 
laid the groundwork for this step in the discussion of Least Squares (Section 5.1). In particular, this Least 
Squares projection equation (5.14) stipulates that 


P=e (e:Te1) ex? and P= ea(e2T ex) ea7 


As e,7e, = e;7 e, = 0 these formulas can not possibly be correct. Returning to the Least Squares discussion 
we realize that it was, perhaps implicitly, assumed that all quantities were real. At root is the notion of the 
length of a complex vector. It is not the square root of the sum of squares of its components but rather 
the square root of the sum of squares of the magnitudes of its components. That is, recalling that the 
magnitude of a complex quantity z is 22, 


(|| ex |)? Aer? e1 rather (|| e1 ||)” 4 ae 


Yes, we have had this discussion before, recall complex numbers, vectors, and matrices (Section 7.1). The 
upshot of all of this is that, when dealing with complex vectors and matrices, one should conjugate before 
every transpose. Matlab (of course) does this automatically, i.e., the ’ symbol conjugates and transposes 
simultaneously. We use x” to denote ‘conjugate transpose’, i.e., 


gaa 
All this suggests that the desired projections are more likely 
Pi =E1 (e:¥e1) ex and P, _ eo(ep#en) en# (9.22) 


Please check that (9.22) indeed jives with (9.21). 
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9.6 The Eigenvalue Problem: Exercises’ 


9.6.1 Exercises 


1. Argue as in Proposition 1 (p. 97) in the discussion of the partial fraction expansion of the transfer 
function that if 7 A k then D,;P, = P;D, = 90. 

2. Argue from this equation (9.17) from the discussion of the Spectral Representation that D;P; = 
P;D; = Dj. 

3. The two previous exercises come in very handy when computing powers of matrices. For example, 
suppose B is 4-by-4, that h = 2 and m, = m2 = 2. Use the spectral representation of B together with 
the first two exercises to arrive at simple formulas for B? and B®. 

4. Compute the spectral representation of the circulant matrix 


6 4 
8 6 
2 8 
4 2 


aon *, WY 
a -, WY CO 


Carefully label all eigenvalues, eigenprojections and eigenvectors. 


6This content is available online at <http://cnx.org/content /m10494/2.3/>. 


Chapter 10 


The Symmetric Eigenvalue Problem 


10.1 The Spectral Representation of a Symmetric Matrix’ 


10.1.1 Introduction 


Our goal is to show that if B is symmetric then 


e each A, is real, 
e each P; is symmetric and 
e each D,; vanishes. 


Let us begin with an example. 


Example 10.1 
The transfer function of 


gee 
B=/111 
i ea 
is 
s-2 1 1 
1 
R(s) 5(8—3) 1 s—2 1 
1 1 8-2 
_{ 2/8 -u8 -1/3 ta 4731/8 
1 
RES) ls 2/3. =1/3) | os | 13. 1/8.4/3 
-1/3 -1/3 -1/3 1/3 1/3 1/3 
1 1 
R(s)= sa Ne 


s—r\ s— 92 


and so indeed each of the bullets holds true. With each of the D; falling by the wayside you may also expect 
that the respective geometric and algebraic multiplicities coincide. 
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10.1.2 The Spectral Representation 


We have amassed anecdotal evidence in support of the claim that each D; in the spectral representation 


h h 
B=) °\jP) +5) D; (10.1) 
j=l j=l 


is the zero matrix when B is symmetric, ie., when B = B’, or, more generally, when B = B” where 
ar : % Od ‘i F 
B# = B’ Matrices for which B = B# are called Hermitian. Of course real symmetric matrices are 


Hermitian. 
Taking the conjugate transpose throughout (10.1) we find 


h h 
BESS GP) +> De (10.2) 
j=l j=l 


That is, the Aj are the eigenvalues of B” with corresponding projections P,;# and _ nilpotents D;" Hence, 
if B = B", we find on equating terms that 


dj = Aj 
R= PM 
and 
D; =D," 


The former states that the eigenvalues of an Hermitian matrix are real. Our main concern however is with 
the consequences of the latter. To wit, notice that for arbitrary 7x, 


(|| Day II) == alt (Diet) pial 
(|| Dees II) =n cD he De 
(|| Dy" *x |)” = 29 DD a 


(|| Dj™~12 ||)” =0 


As DPMS = 0 for every «x it follows (recall this previous exercise (list, item 3, p. 42)) that D;™3-" = 0. 
Continuing in this fashion we find Dat = 0 and so, eventually, D; = 0. If, in addition, B is real then as 
the eigenvalues are real and all the D; vanish, the P; must also be real. We have now established 
Proposition 10.1: 
If B is real and symmetric then 


h 
B=S 7 A;P; (10.3) 
j=l 


where the \,; are real and the P; are real orthogonal projections that sum to the identity and 


whose pairwise products vanish. 
Proof: 
One indication that things are simpler when using the spectral representation is 


h 
Ba Nag By (10.4) 


j=1 
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As this holds for all powers it even holds for power series. As a result, 


h 
ee ) es P; 
j=l 


It is also extremely useful in attempting to solve 
Ba=b 


for x. Replacing B by its spectral representation and b by Jb or, more to the point by iy P;b we 
find 


h h 

YAP e = SOP 

j=l j=l 
Multiplying through by P, gives \, Px = P,bor Pix = am Multiplying through by the subsequent 

Peete Pyb 
P;’s gives Pjx = per Hence, 
h 
L = eet Pix 


(10.5) 
a ee re 0 
We clearly run in to trouble when one of the eigenvalues vanishes. This, of course, is to be expected. 
For a zero eigenvalue indicates a nontrivial null space which signifies dependencies in the columns 
of B and hence the lack of a unique solution to Bz = b. 
Another way in which (10.5) may be viewed is to note that, when B is symmetric, this previous 
equation (9.15) takes the form 


h 
1 
Be ree (10.6) 


Hence, the solution to Bx = b is 


as in (10.5). With (10.6) we have finally reached a point where we can begin to define an inverse 
even for matrices with dependent columns, i.e., with a zero eigenvalue. We simply exclude the 
offending term in (10.6). Supposing that ;, = 0 we define the pseudo-inverse of B to be 


Let us now see whether it is deserving of its name. More precisely, when b € R (B) we would expect 
that x = Btb indeed satisfies Bx = b. Well 


h-1 i h-1 1 h-1 1 h-1 
BES By sPib= Se x BRib = S- sri Pb = So Bp 
j=l “9 jai J j= 


gmk 
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It remains to argue that the latter sum really is b. We know that 
h 
b=) Pjb , bER(B) 
j=l 


The latter informs us that b | N (BF )c As B = B’, we have, in fact, that b L N(B). As Py is 
nothing but orthogonal projection onto N (B) it follows that P;,,b = 0 and so B(Btb) = b, that is, 
x = Btbis asolution to Bx = b. The representation (10.4) is unarguably terse and in fact is often 
written out in terms of individual eigenvectors. Let us see how this is done. Note that if 2 € #(P1) 


then z = P,x and so, 
h 


Ba= BPix = S- Aj Pj Pix = A Pix = yx 


j=l 


ie., is an eigenvector of B associated with A,. Similarly, every (nonzero) vector in Jt (P;) is an 
eigenvector of B associated with A;. 

Next let us demonstrate that each element of Jt (P;) is orthogonal to each element of # (Px) 
when j £k. If € R(P;) and y € KR (P,) then 


gly = (Pjx)" Prey =a' P;Pry =0 


With this we note that if {2j,1,0)2,...,2j,n,} constitutes a basis for 9 (P;) then in fact the union 
of such bases, 
{zjp|(lL<j<h) and (1<p<n;)} 


forms a linearly independent set. Notice now that this set has 


h 
dons 
j=l 


elements. That these dimensions indeed sum to the ambient dimension, n, follows directly from 
the fact that the underlying P; sum to the n-by-n identity matrix. We have just proven 


Proposition 10.2: 
If B is real and symmetric and n-by-n, then B has a set of n linearly independent eigenvectors. 
Proof: 
Getting back to a more concrete version of (10.4) we now assemble matrices from the individual 
bases 

E; = {x51,25,2, aes eat 


and note, once again, that P; = E; Gen Ee and so 


h 
asl 
B=) A;E;(Ej"Ej) E;" 
j=l 


Tunderstand that you may feel a little overwhelmed with this formula. If we work a bit harder we can 
remove the presence of the annoying inverse. What I mean is that it is possible to choose a basis for 
each §t (P;) for which each of the corresponding E; satisfy E;" E; = I As this construction is fairly 
general let us devote a separate section to it (see Gram-Schmidt Orthogonalization (Section 10.2)). 


107 


10.2 Gram-Schmidt Orthogonalization’ 


Suppose that M is an m-dimensional subspace with basis 


{a1,-.-,2m} 


We transform this into an orthonormal basis 


{u, ae # +Im} 
for M via 
1. Set y, = 21 and q, = Tal 
2. yo = £2 minus the projection of x2 onto the line spanned by q;. That is 
-1 
y2 = %2- 41 (an) qi’ &2 = 09. nn’ x2 


Set q2 = Thal and Qe => [q1, 92]. 


3. y3 = £3 minus the projection of x3 onto the plane spanned by q; and q2. That is 


-1 
¥3 = &3 — Q2(Q2" Q2) Qo" x3 = %3 — nn x3 


Set q3 = Tall and Q3 = {q1, g2, 93}. Continue in this fashion through step (m). 


e (m) Ym = Xm minus its projection onto the subspace spanned by the columns of Q,,-1. That is 


m-1 
= 
Yn = Lm — 0, Ana Ona) One tai = be Q9j" Bm 
j=l 
Set dm = 4“, To take a simple example, let us orthogonalize the following basis for R°, 
[Ym ll 
1 1 1 
Ly = 0 tg = i w3 = 1 
0 0 1 


1 m=y1 = 21. 


2. Yo = XQ — ag x = 


3. ys =23-UH' 23 = 


Tr 


We have arrived at 


1 0 0 
n= 0 |aQ= 93 = 0 
0 0 1 


. Once the idea is grasped the actual calculations are best left to a machine. Matlab accomplishes this via 
the orth command. Its implementation is a bit more sophisticated than a blind run through our steps (1) 
through (m). As a result, there is no guarantee that it will return the same basis. For example 
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>x=[1 1 1;0 11 3;00 1]; 
>Q=orth (X) 

Q= 

0.7370 -0.5910 0.3280 
0.5910 0.3280 -0.7370 


0.3280 0.7370 0.5910 


This ambiguity does not bother us, for one orthogonal basis is as good as another. Let us put this into 
practice, via (10.8). 
10.3 The Diagonalization of a Symmetric Matrix® 


By choosing an orthogonal basis { q;,, | 1 << k <n; } for each R(P;) and collecting the basis vectors in 


Q; = ( Oa, Borat. Opa ) 
we find that 
nj 
T 
P=00, =) Gane 
k=1 
As a result, the spectral representation (Section 10.1) takes the form 


B pee Or eras 


h n; T 
= ja Ai gk GRU, ke 
This is the spectral representation in perhaps its most detailed dress. There exists, however, still another 


form! It is a form that you are likely to see in future engineering courses and is achieved by assembling the 
Q; into a single n-by-n orthonormal matrix 


(10.7) 


Q=(a .. @) 


Having orthonormal columns it follows that Q7Q = I. Q being square, it follows in addition that Q7 = Qu. 
Now 

BOj,k = AjQ%,k 
may be encoded in matrix terms via 


BQ=QA (10.8) 


where A is the n-by- n diagonal matrix whose first n, diagonal terms are 1, whose next nz diagonal terms 
are Az, and so on. That is, each A, is repeated according to its multiplicity. Multiplying each side of (10.8), 
from the right, by Q? we arrive at 

B=QAQ™ (10.9) 
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Because one may just as easily write 


QTBQ=A 
one says that Q diagonalizes B. 
Let us return the our example 
111 
B=]/ 11 1 
111 


-1 
a1= 1 
0 
and 
—1 
12> 0 
1 


—1 
2 1 
a 5 
and 
—1 
Te Te 
2 


1 
iB 
1 
and hence 
-V3 -1 V2 
5 1 
Q=(4 qi n )=% v3 -1 v2 
0 oe. age 
and 
0 0 0 
A= 0 0 0 
00 3 
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Chapter 11 


The Matrix Exponential 


11.1 Overview’ 


The matrix exponential is a powerful means for representing the solution to n linear, constant coefficient, 
differential equations. The initial value problem for such a system may be written 


a’ (t) = Az (t) 


x (0) = Xo 


where A is the n-by-n matrix of coefficients. By analogy to the 1-by-1 case we might expect 


a(t) = ety (11.1) 
to hold. Our expectations are granted if we properly define e4* 
each element of At will no suffice? 

There are at least 4 distinct (but of course equivalent) approaches to properly defining e4'. The first 
two are natural analogs of the single variable case while the latter two make use of heavier matrix algebra 
machinery. 


. Do you see why simply exponentiating 


1. The Matrix Exponential as a Limit of Powers (Section 11.2) 

2. The Matrix Exponential as a Sum of Powers (Section 11.3) 

3. The Matrix Exponential via the Laplace Transform (Section 11.4) 

4. The Matrix Exponential via Eigenvalues and Eigenvectors (Section 11.5) 


Please visit each of these modules to see the definition and a number of examples. 

For a concrete application of these methods to a real dynamical system, please visit the Mass-Spring- 
Damper module (Section 11.6). 

Regardless of the approach, the matrix exponential may be shown to obey the 3 lovely properties 

1. £ (e4t) = AeA? = e4tA 

2. eAltitts) — eAti Ate 


3. e4* is nonsingular and (eAt)* = e (AY) 


Let us confirm each of these on the suite of examples used in the submodules. 
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Example 11.1 


If 
1 0 
A — 
0 2 
then 
pat e 0 
0 e2t 
d At et 0 1 0 et 0 
1 dt (e ) => - —" a 
0 2e 0 2 0 
5 efi tte 0 ieee 0 e 60 e? 0 
: 0 e2tit2te J 0 etttetta | 0 eat 0 eta 
-t 
3 (e4t)~" — . 0 — e7 (At) 
0 e @4 
Example 11.2 
If 
0 1 
A = 
-1 0 
then 
“ay cos(t) sin (t) 
e = 
—sin(t) cos (t) 
1. £(e4") = —sin(t) cos (t) sri died = —sin(t) cos (t) 


—cos(t) —sin (t) —cos(t) —sin (¢) 
2. You will recognize this statement as a basic trig identity 


cos(ty +t2) sin (ti + te) cos (ti) sin (€;) cos (tz) sin (t2) 
—sin (t; + t2) cos (ty + te) —sin (ti) cos (ti) —sin (tz) cos (tz) 


3. Using the cofactor expansion we find 


(eat)? _ [ cos(t) —sin(¢) | [ cos(—t) —sin(-t) \ _ oo (At) 
sin(t) cos (t) sin(—t) cos (-—t) 
Example 11.3 
If 
1 
A = 


then 
ee. 
0 1 
0 1 
1. 4 (eAt) = = AeAt 
0 0 
1 ti + te _ 1 ty 1 ta 
0 1 0 1 0 1 
-1 
3. = : = = e— (At) 
0 1 0 il 


11.2 The Matrix Exponential as a Limit of Powers’ 


You may recall from Calculus that for any numbers a and t one may achieve e“ via 
k 
at 
e* = limit (1 + ) 
k—-o0o k 
The natural matrix definition is therefore 
At\* 
e“t = limit (1 + ) 
k-o00 k 


where J is the n-by-n identity matrix. 


Example 11.4 
The easiest case is the diagonal case, e.g., 


for then 
and so (recalling (11.2) above) 


Note that this is NOT the exponential of each element of A. 
Example 11.5 
As a concrete example let us suppose 

0 1 

-1 0 


A- 
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(11.2) 


(11.3) 
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From 
Lie* 6 
I+At= 
—t il 
2 t t t? 
2 xt 1 —t 1 = t 1 t? 
"Di = de 
3 t? e 
(1444 _{ 1-§ t-§& 
3 t? 
: “tog, Ag 
4 —3t? tt 3 
4 — pike —3t? 4 tt Ag 
" 16 8 256 
At\?® sae ae ie sR eee ae 
i+ 7a 44 288 e = | ty 
Mi) 3125 Bh eF25 ae 


We discern a pattern: the diagonal elements are equal even polynomials while the off diagonal 


elements are equal but opposite odd polynomials. The degree of the polynomial will grow with k 
and in the limit we ’recognize’ 


re cos(t) sin (t) 
e — 
—sin(t) cos (t) 
Example 11.6 
If 
1 
A — 
0 0 
then 
At\* ae 
r++) = 
( k 01 
for each value of k and so 
1 ¢ 
eAt — 
01 


11.3 The Matrix Exponential as a Sum of Powers’ 


You may recall from Calculus that for any numbers a and ¢ one may achieve e” via 


love) k 
t=) “ (11.4) 
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The natural matrix definition is therefore 


o°. (At)* 
ee (ar (11.5) 
k=0; Yes 
where A° = I is the n-by-n identity matrix. 
Example 11.7 
The easiest case is the diagonal case, e.g., 
Hie 1 0 
0 2 
for then 
0 
(At)" : 
0 (2t) 
and so (recalling (11.4) above) 
t 
ost e 0 
0 e2t 


Note that this is NOT the exponential of each element of A. 


Example 11.8 
As a second example let us suppose 


and so 
2 3 G 
oat -o4+4-.--: t-R+h-... cos(t) — sin (¢) 


3 5 . 
—t+$-4+... 1-5+4-... —sin (t) cos (t) 
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Example 11.9 


If 
A — 
then 
A? — A? = AK 
and so 
ett =T4+tA= 


11.4 The Matrix Exponential via the Laplace Transform’ 


You may recall from the Laplace Transform module® that may achieve e™ via 


tat / : ) (11.6) 


The natural matrix definition is therefore 
eft = 7} ((sr- 4)"’) (11.7) 


where J is the n-by-n identity matrix. 


Example 11.10 
The easiest case is the diagonal case, e.g., 


1 
A = 
2 
for then 
1 
0 
(sf— A) '=| * ; 
0 s—2 
and so (recalling (11.6) above) 
whe es (4) 0 a. @ 
0 £1 (4) 0 
Example 11.11 
As a second example let us suppose 
0 1 
A — 
-1 0 


and compute, in matlab, 
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> inv(s*eye(2)-A) 


ans = [ s/(s72t+1), 1/(s7*2+1)] 
[-1/(s72+1), s/(s72+1)] 


> ilaplace(ans) 


ans = [ cos(t), sin(t)] 
[-sin(t), cos(t)] 


Example 11.12 
If 


then 


>> inv(s*eye(2)-A) 


ans = [ 1/s, 1/s72] 
[ 0, 1/s] 


> ilaplace(ans) 


, ti 


ans = 1 
Oo, 1] 


11.5 The Matrix Exponential via Eigenvalues and Eigenvectors’ 


In this module we exploit the fact that the matrix exponential of a diagonal matrix is the diagonal matrix 
of element exponentials. In order to exploit it we need to recall that all matrices are almost diagonalizable. 
Let us begin with the clean case: if A is n-by-n and has n distinct eigenvalues, 4;, and therefore n linear 


eigenvectors, s;, then we note that 


As; = \j8; ; j€ {l,...,n} 


®This content is available online at <http://cnx.org/content /m10680/2.13/>. 
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may be written 
A= SNS? (11.8) 
where S = ( By. WSS nseu 8 ) is the full matrix of eigenvectors and A = diag (Ai, A2,..-,An) is the 


diagonal matrix of eigenvalues. One cool reason for writing A as in (11.8) is that 
AP EGS SAS = SAr Se 


and, more generally 

AP = SA*S 
If we now plug this into the definition in The Matrix Exponential as a Sum of Powers (Section 11.3), we 
find 

eft = Sedtg-t 


where e“ is simply 


. Ait Ast Ant 
diag (e°1 6?" po. ce ) 
Let us exercise this on our standard suite of examples. 


Example 11.13 


If 
A = 
then 
S=IA=A 
and so e4! = eAt. This was too easy! 
Example 11.14 
As a second example let us suppose 
0 1 
A — 
-1 0 
and compute, in matlab, 
> [S, Lam] = eig(A) 
S = 0.7071 0.7071 
0 + 0.70711 0 - 0.70711 
Lam = 0 + 1.0000i 0 
0 0 - 1.0000i 
> Si = inv(S) 
Si = 0.7071 0 - 0.70711 


0.7071 0 + 0.70711 


> simple (S*diag(exp(diag(Lam) *t) ) *Si) 


ans = [ cos(t), sin(t)] 
[-sin(t), cos(t)] 


Example 11.15 
If 


then matlab delivers 


> [S, Lam] = eig(A) 


ll 
[o) 
(o) 


Lam 


So zero is a double eigenvalue with but one eigenvector. Hence S is not invertible and we can not 
invoke ((11.8)). The generalization of ((11.8)) is often called the Jordan Canonical Form or the 
Spectral Representation (Section 9.4). The latter reads 


h 
A= S°\jP;+D; 
j=l 
where the 4, are the distinct eigenvalues of A while, in terms of the resolvent R(z) = (zI — Ay, 
4 ie 
J Qri 


R(z) dz 


is the associated eigen-projection and 


1 
D,= — | R(z)(¢-+A;) dz 
= 55 [| ROE) 
is the associated eigen-nilpotent. In each case, C; is a small circle enclosing only Aj. 
Conversely we express the resolvent 


where 
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with this preparation we recall Cauchy’s integral formula (Section 8.2) for a smooth function f 


faa [ Pa 


2771 zZ-a 


where C'(a) is a curve enclosing the point a. The natural matrix analog is 
—1 
f(A) =F [FOR 
Th 
where C'(r) encloses ALL of the eigenvalues of A. For f (z) = e*' we find 
h mj—-1 tk 
At _ vjt [ p. pk 
ett = Se (*+ S- =) (11.9) 


with regard to our example we find, h = 1, A; = 0, P) =I, m,; = 2, D; = Aso 


Let us consider a slightly bigger example, if 


then 


> R = inv(s*eye(3)-A) 


R=[ 1/(s-1), 1/(s-1)72, 0] 
[ 0, 1/(s-1), 0] 
[ 0, 0, 1/(s-2)] 


and so A, = 1 and A» = 2 while 


and so m, = 2 


and 


Py 


II 
x a a =) 
e- Oo O&O 
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and m2 = 1 and Dj = 0. Hence 


et = et (P; +tD,) + e7' P, 


e tet 
0 ef 
0 O. eé 


11.6 The Mass-Spring-Damper System’ 


ay k 2c 


Figure 11.1: Mass, spring, damper system 


If one provides an initial displacement, x9, and velocity, vp, to the mass depicted in Figure 11.1 then one 
finds that its displacement, x (t) at time t satisfies 


Po) | dex (t) 


qe + 2¢—y + ka (t) = 0 (11.10) 


x (0) = 2x0 


x’ (0) = V0 


where prime denotes differentiation with respect to time. It is customary to write this single second order 
equation as a pair of first order equations. More precisely, we set. 


and note that (11.10) becomes 
muy’ (t) = (— (kui (t))) — 2cuz (t) (11.11) 
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ur’ (t) = ue (t) 


T 
Denoting u (t) = ( uz (t) up (t) ) we write (11.11) as 


i 
uw(t)=4ut) , A=[ (11.12) 


We recall from The Matrix Exponential module that 


u(t) = e4*u (0) 


We shall proceed to compute the matrix exponential along the lines of The Matrix Exponential via Eigen- 
values and Eignevectors module (Section 11.5). To begin we record the resolvent 


Oy -1 2c+mz m 
1 mz + 2ez +k _k hie 


The eigenvalues are the roots of mz? + 2cz +k, namely 


eid 
Ses d= Vc? —mk 


dy = zeae , d= Vc?—mk 
m 


We naturally consider two cases, the first being 


1. d#0. In this case the partial fraction expansion of R(z) yields 


Pon -1 1 d-—c —m ze -1 1 c+d om = 1 py “1 p 
poe Od ie: Gated z— 2X» 2d age Beep BM Yo eg 
and so e4! = e*1' P, + €2' Py. If we now suppose a negligible initial velocity, i.e., vo, it follows that 
a(t) = st (e* (d—0) + &* (c+-d)) (11.13) 


If d is real, ie., if c? > mk, then both \; and Xz are negative real numbers and z(t) decays to 0 
without oscillation. If, on the contrary, d is imaginary, i.e., c? < mk, then 


a(t) =e") (cos (\d|t) + qe (at) (11.14) 


and so x decays to 0 in an oscillatory fashion. When (11.13) holds the system is said to be overdamped 
while when (11.14) governs then we speak of the system as underdamped. It remains to discuss the 
case of critical damping. 


. d= 0. In this case \y = Ag = —4/ -. and so we need only compute P; and D,. As there is but one 


P; and the P; are known to sum to the identity it follows that P,; = I. Similarly, this equation (9.16) 
dictates that 


Di, =AP,-yP=A-AI= 


| 
se 3 
a 


On substitution of this into this equation (11.9) we find 
1+t,/£ t 
@ /* ) 1—t,/* 


Under the assumption, as above, that vp = 0, we deduce from (11.15) that 


x(t) = oe (tVe) (: 3 i) Xo 
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Chapter 12 


Singular Value Decomposition 


12.1 The Singular Value Decomposition’ 


The singular value decomposition is another name for the spectral representation of a rectangular matrix. 
Of course if A is m-by-n and m 4 n then it does not make sense to speak of the eigenvalues of A. We 
may, however, rely on the previous section to give us relevant spectral representations of the two symmetric 
matrices 


e ATA 
e AAT 


That these two matrices together may indeed tell us ’everything’ about A can be gleaned from 


WN (ATA) = W (A) 
W (AAT) =. (AT) 
R (ATA) =R(AT) 


R (AAT) = R(A) 
You have proven the first of these in a previous exercise. The proof of the second is identical. The row and 
column space results follow from the first two via orthogonality. 
On the spectral side, we shall now see that the eigenvalues of AA’ and A’ A are nonnegative and that 
their nonzero eigenvalues coincide. Let us first confirm this on the adjacency matrix associated with the 
unstable swing (see figure in another module (Figure 3.2: A simple swing)) 


0 10 0 
A=] -1 0 1 0 (12.1) 
Oe.) OT 


The respective products are 


100 
AAT=] 0 2 0 
001 
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ATA= 


| 
ra 
oor SO 
oro 
FOO Oo 


Analysis of the first is particularly simple. Its null space is clearly just the zero vector while 4; = 2 and 
A2 = 1 are its eigenvalues. Their geometric multiplicities are ny = 1 and ng = 2. In A? A we recognize 
the S matrix from the exercise in another module?and recall that its eigenvalues are X41 = 2, A2 = 1, and 
3 = 0 with multiplicities ny = 1, no = 2, and ng = 1. Hence, at least for this A, the eigenvalues of AA? 
and A’ A are nonnegative and their nonzero eigenvalues coincide. In addition, the geometric multiplicities 
of the nonzero eigenvalues sum to 3, the rank of A. 


Proposition 12.1: 

The eigenvalues of AA? and A™ A are nonnegative. Their nonzero eigenvalues, including geometric 
multiplicities, coincide. The geometric multiplicities of the nonzero eigenvalues sum to the rank of 
A. 

Proof: 

If AT Ax = Ax then x7 AT Ax = dx™ 2, ice., (|| Ax ||)” = (|| x ||)? and so \ > 0. A similar argument 
works for AA’. 


Now suppose that \; > 0 and that {2;,,};%, constitutes an orthogonal basis for the eigenspace R (Pj). 


Starting from 
A? Ax; = AjXjk (12.2) 


we find, on multiplying through (from the left) by A that 
AA? Azjx = Aj AX;,b 


ie., A; is an eigenvalue of AA? with eigenvector Ax;,, so long as Axj,, 4 0. It follows from the first 
paragraph of this proof that || Av; , ||= ./;, which, by hypothesis, is nonzero. Hence, 


l<k<n, (12.3) 


is a collection of unit eigenvectors of AA? associated with \;. Let us now show that these vectors are 
orthonormal for fixed 7. 
Ys Vink = stheAT Aas = UF tik =0 
We have now demonstrated that if \; > 0 is an eigenvalue of A? A of geometric multiplicity n;. Reversing the 
argument, i.e., generating eigenvectors of A’.A from those of AA? we find that the geometric multiplicities 
must indeed coincide. 

Regarding the rank statement, we discern from (12.2) that if \; > 0 then aj, € R (A7A). The union of 
these vectors indeed constitutes a basis for R (ATA), for anything orthogonal to each of these x;,, necessarily 
lies in the eigenspace corresponding to a zero eigenvalue, i.e., in VY (ATA). AsR (ATA) =R (AT) it follows 
that dimR (AA) =r and hence the nj, for A; > 0, sum to r. 

Let us now gather together some of the separate pieces of the proof. For starters, we order the eigenvalues 
of A? A from high to low, 

Ay > Ag >+++>An 
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and write 
AP A=XA,XT (12.4) 


where 
X= {X1,..., Xn} 5 X; = (ipiyicig ting 


and A, is the n-by-n diagonal matrix with \, in the first n, slots, Ag in the next ng slots, etc. Similarly 
AA SR (12.5) 


where 
Y= Lig ceca Ya ’ Y; = 1U5 peng Binet 


and Aj, is the m-by-m diagonal matrix with A; in the first n; slots, Ag in the next mg slots, etc. The yj; x 
were defined in (12.3) under the assumption that A; > 0. If A; = 0 let Y; denote an orthonormal basis for 


N (AAT). Finally, call 
o=V% 


and let © denote the m-by-n matrix diagonal matrix with o, in the first n, slots, o2 in the next ng slots, 
etc. Notice that 
od = An (12.6) 
yt = Ay (12.7) 
Now recognize that (12.3) may be written 
AX; k = OFU5,k 
and that this is simply the column by column rendition of 


AX =Y™% 


As XX? =I we may multiply through (from the right) by X7 and arrive at the singular value decom- 
position of A, 
A=YuXT (12.8) 


Let us confirm this on the A matrix in (12.1). We have 


oO FF Oo O&O 


oO fF Oo - 
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Hence, 
v2 0 0 0 
n= 0 100 (12.9) 
0 0 1 0 
and so A= YX" says that A should coincide with 
= 9 Lt Qg 
0 1.0 v2 0 0 0 va ve 
0 1 0 O 
1 0 0 0 1 0 0 
0 0 0 1 
00 1 0 0 1 0 : ; 
oe 0 a 0 


This indeed agrees with A. It also agrees (up to sign changes on the columns of X) with what one receives 
upon typing [Y, SIG, X] = scd(A) in Matlab. 
You now ask what we get for our troubles. I express the first dividend as a proposition that looks to me 
like a quantitative version of the fundamental theorem of linear algebra. 
Proposition 12.2: 
If YXX7 is the singular value decomposition of A then 


1. The rank of A, call it r, is the number of nonzero elements in &. 

2. The first r columns of X constitute an orthonormal basis for R (A’). The n—rlast columns 
of X constitute an orthonormal basis for -V (A). 

3. The first r columns of Y constitute an orthonormal basis for R (A). The m—r last columns 
of Y constitute an orthonormal basis for (A?) 


Let us now ’solve’ Ax = b with the help of the pseudo-inverse of A. You know the ’right’ thing to do, 
namely reciprocate all of the nonzero singular values. Because m is not necessarily n we must also be careful 
with dimensions. To be precise, let U* denote the n-by-m matrix whose first n, diagonal elements are oe 


whose next nz diagonal elements are - and so on. In the case that o, = 0, set the final n, diagonal elements 
of 4+ to zero. Now, one defines the pseudo-inverse of A to be 


At=xyty™ 
In the case of that A is that appearing in (12.1) we find 


V2 0 0 
sea) © 1.8 
0 O11 
0 OO 
and so 
4 1 1 
mm Oi ae 0 a A af 
0 1 0 0 0 1 0 
At = 1 0 0 
0 Oo 0 1 0 Oot 
i ‘ 001 
we Oa 0 O 0 
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therefore, 


BR 


At= 


oO oO FF CO 
On Oo || 
te <> ee => a 


in agreement with what appears from pinv(A). Let us now investigate the sense in which A?* is the inverse 
of A. Suppose that b € R™ and that we wish to solve Ar = b. We suspect that Atb should be a good 
candidate. Observe by (12.4) that 


(ATA) ATb=XA,XTXETYTD 
because X7 X = I 
(A? A) Atb= XA, ZUtY7D 


by (12.6) and (12.7) 
(APA) AtTb=XETIEUTYTS 


because DP UH+ = NF 
(A™A) Atb= X=TY7) 


by (12.8 
mre (ATA) Ath = ATD 


that is At satisfies the least-squares problem A? Ax = A” b. 
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Glossary 


A. A Basis for the Column Space 


Suppose A is m-by-n. If columns {c; | 7 =1,...,r} are the pivot columns of A;eq then columns 
{c; | j =1,...,r} of A constitute a basis for Ra (A). 


A Basis for the Left Null Space 


Suppose that A? is n-by-m with pivot indices {c; | j = {1,...,r}} and free indices 
{ce; |j ={r+1,...,m}}. A basis for ” (A) may be constructed of m — r vectors 
{z1,2?,...,2™-"} where z*, and only z”, possesses a nonzero in its c,4, component. 


A Basis for the Null Space 


Suppose that A is m-by-n with pivot indices {c; | j = {1,...,r}} and free indices 
{cj |j={r4+l1,...,n}}. A basis for (A) may be constructed of n — r vectors 
{z1,2?,...,2"-"} where z*, and only z*, possesses a nonzero in its c,+, component. 


A Basis for the Row Space 
Suppose A is m-by-n. The pivot rows of Ayeq constitute a basis for RA (AT). 
A subset S of a vector space V is a subspace of V when 
1. if « and y belong to S then so does x + y. 
2. if x belongs to S and t is real then ta belong to S. 


axial resistance 


B Basis 


Any linearly independent spanning set of a subspace S$ is called a basis of S. 


C Capacitance of a Single Compartment 


C= nase 
capacitance of cell body 
Cob = Acne 
Column Space 


The column space of the m-by-n matrix S is simply the span of the its columns, i.e. 
Ra(S)={Sa|ae€ R"}. This is a subspace (Section 4.7.2) of R”. The notation Ra stands for 
range in this context. 


Complex Addition 
Z1 + Z2 —a ah + x2 +4 (yi + y2) 
Complex Conjugation 


A= 21 — ty, 
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Complex Division 


1 — 21 22 __ Ui tet+yryo2ti(eey1—21y2) 
22” «22 BQ x27 +y2 


Complex Multiplication 


2122 = (@1 + iy) (fo + ty2) = 1122 — Yryo +7 (@1y2 + Loy1) 


D Dimension 


The dimension of a subspace is the number of elements in its basis. 


E Elementary Row Operations 
1. You may swap any two rows. 


2. You may add to a row a constant multiple of another row. 


F force balance 


1. Equilibrium is synonymous with the fact that the net force acting on each mass must vanish. 
2. In symbols, 
yi — y2— fi =0 


Y2— 43 — fz =0 


y3 — ya — fg = 0 


3. or, in matrix terms, By = f where 


a i ra 8 
f= |) and B=]0 1 —-1 O 
fa Ge Gri Het 


H Hooke’s Law 


1. The restoring force in a spring is proportional to its elongation. We call the constant of 
proportionality the stiffness, k;, of the spring, and denote the restoring force by y;. 
2. The mathematical expression of this statement is: y; = kje,;, or, 


3. in matrix terms: y = Ke where 


ky 0 O O 

0 kz, O O 
K= 

0 O kz O 

0 O O kg 


I Inverse of S 


Also dubbed "§ inverse" for short, the value of this matrix stems from watching what happens 
when it is applied to each side of Sx = f. Namely, 


(SSS Ses Saas See) 
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Invertible, or Nonsingular Matrices 
Matrices that do have an inverse. 


Example: The matrix S that we just studied is invertible. Another simple example is 


L Left Null Space 


The Left Null Space of a matrix is the null space (Section 4.2) of its transpose, i.e., 
MN (AT) ={yeR™ | ATy=0} 


Linear Independence 


A finite collection {s1, 82,...,8n} of vectors is said to be linearly independent when the only 
reals, {1,22,...,2n} for which v1 + v2 +--+: +2, =O are v1 = 4 =--- = 2, = 0. In other 
words, when the null space (Section 4.2) of the matrix whose columns are {81, 52,...,5n} 
contains only the zero vector. 


M Magnitude of a Complex Number 


|za| = Jaa = far? + 1? 


membrane resistance 


Roo 
ray 
N Nernst potentials 
RT, [Nal RT, [kK] 
Ena = —1 ° d Ex = —1 ° 
een Males TR: 
where R is the gas constant, T’ is temperature, and F' is the Faraday constant. 


Null Space 
The null space of an m-by-n matrix A is the collection of those vectors in R” that A maps to the 


zero vector in R™. More precisely, 


WN (A) ={a €R"| At =0} 


O orthogonal projection 


A matrix P that satisfies P? = P is called a projection. A symmetric projection is called an 
orthogonal projection. 


P Pivot Column 
Each column of A;.q that contains a pivot is called a pivot column. 
Pivot Row 
Each nonzero row of Ayeq is called a pivot row. 


Pivot 
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The first nonzero term in each row of A;eq is called a pivot. 
poles 


Also called singularities, these are the points s at which Ly, (s) blows up. 


R Rank 
The number of pivots in a matrix is called the rank of that matrix. 
resistance of cell body 
Rep = AcbPm- 
Row Space 


The row space of the m-by-n matrix A is simply the span of its rows, i.e., 


a(A?) ={ATy|yeR™} 


S singular matrix 
A matrix that does not have an inverse. 


Example: A simple example is: 


Span 


A finite collection {s1,s2,...,5n} of vectors in the subspace S' is said to span S if each element of 
S can be written as a linear combination of these vectors. That is, if for each s € S' there exist n 
reals {x1,%2,...,%y} such that s = 7181 + vaso +--+ + an Sn. 


T The Row Reduced Form 


Given the matrix A we apply elementary row operations until each nonzero below the diagonal is 
eliminated. We refer to the resulting matrix as Ajeq. 


V_ Voltage-current law obeyed by a capacitor 


The current through a capacitor is proportional to the time rate of change of the potential across 
it. 
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